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I. INTRODUCTION 



The 1973 oil crisis, where the price of a barrel of oil 
jumped higher than 210 percent and affected every aspect of 
the world economy , direct ly affected the shipping transporta- 
tion industries . Whereas the price of oil was formerly of 
modest importance, it became a prime concern . Extensive asso- 
ciation to the emerged problem has to do with the design of 
autopilots for ships, which can minimize propulsive losses 
caused by added resistance due to steering. 

One motivation for the design of an optimum steering 
controller was the work done by Nomoto and Motoyama who 
claimed that reduction in propulsion loss could amount to a 
1 percent savings in fuel consumption. 

For most commercial ships a 1 to 2 percent saving in 
fuel costs , justifies the expense of fitting an autopilot 
which has the capability of producing this savings [Ref. 1] . 

Chapter 2' addresses what type of computer model can be 
used to represent the ship. Several of these models are 
investigated , such as simulation from the equations of motion 
and the Nomoto third order which was developed from the 
equation of motion. 

Chapter 3 addresses the problem of selecting an adequate 
cost function to represent the added resistance due to 
steering . 

The problem of finding the best controller design to 
provide a minimum value of added resistance due to steering, 
for regular and irregular seas, is studied in Chapters 4 and 
5 respectively by using as a model of the ship the equations 
of motion and a function minimization subroutine. 

Chapter 6 indicates how the controller parameters can be 
adjusted in any encounter environmental condition by means 
of an adaptive control. 
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Chapter 7 presents the conclusions from the experimental 
work . 
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II. DESCRIPTION OF COMPUTER MODELS 



The most accurate model which can represent the ship/ 
steering dynamics is the model which is based upon the equa- 
tions of motion as defined by series expansion including all 
terms (both linear and nonlinear) . 

Using experimentally measured hydrodynamic coefficients, 
for the SL-7 high speed containership , a computer program 
was developed in order to provide a computer simulation for 
the ship. Figure 2.1 shows the block diagram and the computer 
program can be found in Appendix A [Ref. 2,3]. 




Figure 2.1 Block Diagram of Ship and Control System 
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Figure 2.2 Determination of Nomoto Third Order Model 

As a next choice, from the equations of motion we derive 
the Nomoto third order transfer function. Figure 2.2 shows 
the block diagram. 

Using the scheme of Figure 2.2 which includes the func- 
tion minimization subroutine with both yaw and sway equa- 
tions, with the linear terms only we obtain the appropriate 
coefficients of the Nomoto third model. (Including in the 
equations of motion nonlinear terms we can see that the 
perturbations were small enough). 

The function minimization subroutine used was BOXPLX, 
which was programed by R. Hilleary.The task of the already 
mentioned subroutine is to find the minimum of any function 
and is subjected to explicit constraints of the variables or 
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implicit 
tion it 
The 

got by 
for 16, 



constraints on functions of the variables. In addi- 
can handle a maximum of 25 variables [Ref. 4]. 
extracted results were very close to the results we 
the analytic solution and are tabulated below only 
23 and 32 knots. 
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Nomoto Third Order Model 
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III. STEERING PERFORMANCE CRITERION 



Studying the literature, one can see many approaches to 
the problem of optimizing an automatic ship steering 
controller for maximum reduction in fuel consumption. Since 
added resistance is directly related to both rudder activity 
and yawing motion we can express a measure of this added 
resistance given as a performance criterion with the 
formula : 



• Lambda(^ )=weighting factor 

This formula defines an approximate drag due to steering 
for small amplitude oscillations about a steady-state pivot 
point of the ship during yawing at the natural frequency of 
the closed-loop ship/ steering system . (About 0.05 rad/sec for 
the SL-7 containership ) . It is also convenient for shipboard 
use because yaw error and rudder angle can be easily meas- 
ured . 

During this study the values for the weighting factors 
are taken from R. E. Reid's work for the high-speed contain- 
ership SL-7 [Ref. 5]. 




(3.1) 



Where : 

• Delta (S )=rudder angle 
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TABLE 2 

Weighting Factors 



Ship's Speed, Knots Lambda 

16 16.796 
23 8.128 
32 4.2 



In Table 2 weighting factors for the operating range of 
the ship are tabulated. 
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IV. CONTROLLER DESIGN FOR REGULAR SEAS 



Our goal now, is to estimate the system's performance in 
a seaway. To accomplish the above task first we must deter- 
mine a suitable representation of the external disturbances 



ciently accurate computer ship model and a steering perform- 
ance criterion have been defined. 

It can be postulated that a sufficiently accurate model 
of the seaway itself, is a representative modeling of forces 
and moments exerted on the ship [Ref. 6,7,8]. An explana- 
tion of what a sea comprises, and how predicted or observed 
sea states can be analyzed to determine the forces and 
motions of a body in a sea was studied by Michael [Ref. 9]. 

In this chapter we will use as a seaway representation 
the regular sea model in which the forces on the ship 
[Ref. 10], have the form: 



Values for the exciting force ( ) for different 

encounter angles and different encounter frequencies as 



spondence between wave height and sea state which were used 
in this work. 



on the ship by the sea. This can be done, since a suffi- 




(4.1) 



Where : 

• Re =exciting force 

• CJe ^encounter frequency 

• w h =wave height 




i =phase angle 



well, were taken from [Ref. 3]. Table 3 shows the corre- 
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The controller used in the present study is shown in 
Appendix B and its form is repeated here in Figure 4.1 . 



TABLE 3 

Sea State Vs Wave Height 



Sea State 

(Beaufort scale) 
1 
2 
4 
6 
7 



Wave Height 

(in feet) 
4.0 

4.5 

5.5 
10.0 
17.5 




Figure 4.1 Controller C 

Values for the optimal gains for the above controller 
and values for the cost J as well, are shown below in Tables 
4,5,6,7,8,9,10,11 

These values were obtained using: 

• Sea states 1-2-4-6-7 (Beaufort scale). 

• Encounter angles 0 0 - 30° - 60 0 - 90 0 . 

• Encounter frequencies 0.2-0.4-0.6-0.75-1.5 rad per sec. 

• Constant speed 23 knots. 

A careful analysis of the extracted results lead us to 
conclude : 
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• The maximum deviation of controller parameter values 
occured at 0.4 rad per sec encounter frequency for all 
used encounter angles and sea states. 

• For the same encounter f requency , f or all encounter 
angles studied, the controller parameter values 
increasing smoothly as the sea state increases. 

• For all encounter frequencies and for specific 

encounter angles, we observe that the cost increases at 
higher sea states. 

• For 0.4 rad per sec encounter frequency the cost 
changes rapidly for sea state 6 and 7 as we move 
through the encounter angles . 

• For the same encounter angles and sea states the cost 
decreases at higher encounter frequencies. 

• For all encounter angles and encounter frequencies the 
maximum cost occured for sea state 7. 

• The cost had it's maximum value for encounter angle 90° 
and encounter frequency 0.4 rad per sec. 

Furthermore in order to obtain the behavior of the 
rudder and yaw motion of the ship, transient response plots 
were obtained for controller 'C' at ship's speed 23 knots, 
different sea states and encounter angles as shown in 
Figures 4.2 through 4.15 

These plots were obtained using the program of Appendix 
C. 

From these plots it is verified that for the same 
encounter frequencies and encounter angles , rudder and yaw 
motion increases in amplitude as the sea states rises. 

Finally Figures 4.16,4.17,4.18,4.19,4.20, show the 
results of an experiment in which we were changing one 
controller parameter keeping the rest fixed, and we observed 
that the cost does not change significantly in the vicinity 
of the actual value. 
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Using that as a fact we can postulate that for a 
specific sailing mode for the controller parameter values 
high accuracy is not required. 

Comparing controller C with controller A it is obvious 
that the parameter values of controller A vary over wider 
range . 
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TABLE 4 

Controller C For 0° Encounter Wave Angle 



Sea State: 
Kl = 

Tl= 

T2 = 

T3 = 

T4 = 

Cost J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

Cost J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

Cost J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

Cost J= 



Sea State: 
Kl= 

Tl = 

T2 = 

T3 = 

T4= 

Cost J = 



Encounter Frequency 0.2 Rads 



Per Sec 



1 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

0.6382236E-33 

Encounter Frequency 



2 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 8077518E-33 

0.4 Rads Per Sec 



1 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

O.9306967E-35 

Encounter Frequency 



2 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 1177913E-34 

0.6 Rads Per Sec 



1 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 1868590E- 36 



2 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 2364934E-36 



Encounter Frequency 0.75 Rads Per Sec 



1 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0. 1917420E-38 

Encounter Frequency 



2 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 2426735E- 38 

1.5 Rads Per Sec 



1 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0. 1476825E-38 



2 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0. 1869107E-38 
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TABLE 5 

Controller C For 30° Encounter Wave Angle 



Encounter Frequency 0.2 Rads Per Sec 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

Cost J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

Cost J= 



Sea State: 
Kl= 

Tl= 

T2= 

T3 = 

T4 = 

Cost J= 



Sea State: 
Kl = 

Tl = 

T2 = 

T3 = 

T4 = 

Cost J= 



1.7432022 
35 . 1575165 
22.0137024 
13.5792084 
22.0728760 
0 . 594662E+00 



1.7434111 
35.0142517 
22.0644684 
13.7818375 
22.1337128 
0 . 7502608E+00 



Encounter Frequency 0.4 Rads Per Sec 



0.2120571 
92.9998932 
24.2121429 
39 . 9739990 
32.8156433 
0 . 230148 IE-02 



0.1997837 
99.9923096 
24.2589569 
39.9999542 
34.0418396 
0 . 2901770E- 02 



Encounter Frequency 0.6 Rads Per Sec 



1.8069019 
1.1657065 
12.4824829 
13 . 7605896 
9.2525921 
0 . 1797669E-02 



1.7886524 
1.1832294 
13.3286743 
15.5890350 
9.2333460 
0 . 2275 77 IE- 02 



Encounter Frequency 0.75 Rads Per Sec. 



2.2712851 
0. 84311'* 1 
24.1657867 
22.2913666 
9.8764191 
0. 1094168E-02 



2.2712851 
0.8431141 
24.1657867 
22.2913666 
9.8764191 
0 . 1384705E-02 



Encounter Frequency 1.5 Rads Per Sec 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

Cost J= 



1.6945763 
0.2571034 
25 . 3036194 
25.4920502 
26.4106903 
0 . 5423090E- 04 



1.6945763 
0.2571034 
25.3036194 
25.4920502 
26.4106903 
0 . 6863636E-04 
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TABLE 6 

Controller C For 60° Encounter Wave Angle 



Encounter Frequency 0.2 Rads Per Sec 



Sea State 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

Cost J= 



Sea State: 
Kl = 

Tl= 

T2 = 

T3 = 

T4 = 

Cost J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T4= 

Cost J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

Cost J= 



2 

29 

18 

10 

18 



1659079 

1530609 

6070404 

7567606 

3101044 



0. 7797163E=00 



2.1752367 
28.9752655 
18.4342346 
10.5662384 
18.4618378 
0 . 9817770E+00 



Encounter Frequency 0.4 Rads Per Sec 



0.8408279 
73.1355896 
53.1565094 
22.6398926 
49 . 9999390 

0. 7811032E+000. 978357 1E+00 



0.8623953 
74.4368439 
57.8289948 
24.3586884 
49 . 9900360 



Encounter Frequency 0.6 Rads Per Sec 



1.0592918 
0 . 1000010 
19 . 7355957 
49 . 9999390 
0.2017 35 IE- 02 



0.9676542 
0 . 1343303 
25 . 3748169 
33.4176483 
0. 2553086E-02 



Encounter Frequency 0.75 Rads Per Sec 



2.2673655 
0.8431506 
14.1734161 
11.3151855 
15 . 1259995 
0 . 2704291E-02 



2.2673655 
0.8431506 
14.1734161 
11.3151855 
15 . 1259995 
3422481E-02 



Encounter Frequency 1.5 Rads Per Sec 



Sea State: 
Kl = 

Tl= 

T2 = 

T3 = 

T4= 

Cost J= 



2.3734770 
0.5810375 
8.2907715 
6.7814789 
7.9672575 
0 . 1145959E-02 



2 

2.3808193 
0.5764611 
0220699 
6438484 
3380499 



7 

5 

9 



0. 1450353E-02 
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TABLE 7 

Controller C For 90° Encounter Wave Angle 





Encounter Frequency 


CM 

O 


Rads 


Per Sec 


Sea State: 


1 






2 


Kl = 


1.3473454 






1.3487740 


Tl = 


35.3837585 






35.2085876 


T2 = 


19.5286560 






19.5352631 


T3 = 


15.7849731 






15.8424835 


T4= 


19.3971252 






19.5299683 


Cost J= 


0. 1756000E+01 






0 . 2201721E+01 




Encounter Frequency 


0.4 


Rads 


Per Sec 


Sea State: 


1 






2 


Kl= 


1.6298742 






1.6293240 


Tl= 


31.0220642 






31.2911530 


T2 = 


18.8860321 






19.0679169 


T3 = 


12.3363419 






12.2226353 


T4= 


18.8753967 






19.0280762 


Cost J= 


0 . 6726980E + 00 






0 . 8399093E + 00 




Encounter Frequency 


0.6 


Rads 


Per Sec 


Sea State: 


1 






2 


Kl = 


1.9138527 






1.9098577 


Tl= 


1.1895151 






1.1946430 


T2 = 


10.0705585 






10.1293569 


T3 = 


11.7410755 






11.8040771 


T4= 


10.0737257 






10.1994019 


Cost J= 


0 . 5086016E-01 






0 . 6385970E-01 




Encounter Frequency 


0.75 


» Rads Per Sec 


Sea State: 


1 






2 


Kl= 


1.8993120 






1.8893242 


Tl = 


0.4503812 






0.4850104 


T2 = 


27.3063202 






31.5392303 


T3 = 


21.9325409 






27.4768982 


T4= 


28.9246674 






26.0999603 


Cost J= 


0 . 2837797E-01 






0 . 3586013E-01 




Encounter Frequency 


1.5 


Rads 


Per Sec 


Sea State 










Kl= 


2.4287634 






2.4257565 


Tl= 


0.6342447 






0.6373515 


T2 = 


7.4652443 






7.9898472 


T3 = 


5.9140081 






6.1980715 


T4 = 


5.8591156 






5.1957550 


Cost J= 


0. 1057035E-01 






0. 1336560E-01 
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TABLE 8 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

COST J= 



Sea State: 
Kl = 

Tl = 

T2 = 

T3 = 

T4 = 

COST J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

COST J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

COST J= 



Sea State: 
Kl= 

Tl = 

T2 = 

T3 = 

T4 = 

COST J= 



Controller C For 0° Encounter Wave Angle 



Encounter Frequency 0.2 Rads 



4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 9972246E- 33 



4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0. 1454213E-34 



4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 2919672E-36 



4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 2995968E- 38 



6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 3988896E- 32 



6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 5816855E- 34 



6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 1167868E-35 



6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0. 1198388E-37 



Per Sec 



7 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0. 1221600E-31 

Per Sec 



7 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0. 1781412E-33 

Per Sec 



7 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 3576599E-35 

Per Sec 



7 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 3670063E- 37 



Encounter Frequency 0.4 Rads 



Encounter Frequency 0.6 Rads 



Encounter Frequency 0.75 Rads 



Encounter Frequency 1.5 Rads Per Sec 



4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 2307540E- 38 



1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 9230160E- 38 



7 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 
0 . 2826737E-37 
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TABLE 9 

Controller C For 30° Encounter Wave Angle 



Encounter Frequency 0.2 Rads Per Sec 



Sea State: 
Kl = 

Tl= 

T2 = 

T3 = 

T4= 

COST J= 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

COST J= 



1.7466650 
35.2978973 
22.2485657 
13.6511078 
22.1172638 
0 . 9229971E+00 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

COST J= 



1.7746067 
36.3333435 
23.2952423 
13.4007721 
23.5081635 
0. 3496342E+01 



1.8194151 
48 . 3066101 
45.8587494 
25.3368378 
23 . 1087189 
0 . 9325538E + 01 



Encounter Frequency 0.4 Rads Per Sec 



0.2102537 
92.8943787 
23.9096680 
39 . 9628601 
33 . 1146698 
0 . 3595694E-02 



0.2002187 
99.9424896 
24.3453674 
39.9452667 
33.9409027 
0. 1431070E-01 



0.2060161 
95 . 9424133 
24.2843170 
38.9973450 
33.4644012 
0 . 43769 12E- 0 1 



Encounter Frequency 0.6 Rads Per Sec 



1.8260603 
1.1601276 
11.2977753 
12.5801849 
10.2739115 
0 . 2808505E- 02 



1.8175459 
1.1572495 
10.3958378 
11.7078133 
10.7313347 
0. 1122463E-01 



1.8287125 
1.1652012 
10.5659571 
11.7124157 
10.5683947 
0 . 3429835E-01 



Encounter Frequency 0.75 Rads Per Sec 



Sea State 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

COST J= 



4 

2.2723074 
0.8993683 
22.2714996 
20.4913177 
10.0878286 
0. 1709419E-02 



2.3248615 
0.8716087 
13.4890442 
10.3330994 
12.5009842 
0. 6832119E-02 



2.3073282 
0.8750648 
13.0364380 
10.1076660 
13 . 1554108 
0 . 2090112E- 01 



Encounter Frequency 1.5 Rads Per Sec 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

COST J= 



1.9173651 
0.6582609 
38 . 7910156 
19.3868103 
17 . 9940795 
0 . 8474453E- 04 



1.8726044 

0.3396787 

31.8340454 

29.9999695 

24.4667816 

0.3389677E-03 



1.6678009 
0.3818831 
27 . 1428680 
26.3686218 
27 . 3681793 
0. 1038366E-02 
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TABLE 10 

Controller C For 60° Encounter Wave Angle 



Encounter Frequency 0.2 Rads Per Sec 



Sea State 
Kl= 

Tl = 

T2 = 

T3 = 

T4 = 

COST J= 



Sea State 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

COST J= 



2.1778851 
29 .6119080 
18.7636871 
10.2442350 
18 . 3449402 
0 . 1205 183E+01 



2.2266541 
31.9857025 
19.4135590 
7.8426695 
19.8131866 
0 . 4414552E+01 



2.2981157 
35.2284546 
24.2031403 
8.8281393 
24.7627716 
0. 1080894E+02 



Encounter Frequency 0.4 Rads Per Sec 



0.8723865 
77 . 6901245 
57.4288635 
21.6142426 
58 . 1057587 
0. 1194492E+01 



1.3212109 
1.1219873 
14.2327118 
35.2742767 
14.2909546 
0 . 3620569E+01 



1.7475338 
0.9073439 
18.7083435 
29 . 6310730 
20.4921570 
0. 1005195E+02 



Encounter Frequency 0 . 6 Rads Per Sec 



Sea State: 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

COST J= 



Sea State 
Kl = 

Tl= 

T2 = 

T3 = 

T4 = 

COST J= 



0.9356804 
0.2035089 
25.2695007 
18 . 3522491 
39 . 1636200 
0. 3151959E-02 



1.0065041 
0.1434612 
28 . 2131348 
16.1927643 
29.4095917 
0. 1260864E-01 



0.9630518 
0.1086842 
24.0041809 
17.7048645 
38 . 1658936 
0 . 3862068E- 01 



Encounter Frequency 0.75 Rads Per Sec 



2.2310591 

0.8066239 

18.4929504 

16.4596100 

12.6112747 

0.4225172E-02 



2.2579184 
0.7918448 
14.8840103 
12.1710510 
14.9709930 
0. 1688830E-01 



2.2688522 
0.7898693 
14.6062164 
11.3906097 
14.7650146 
0 . 5162057E-01 



Encounter Frequency 1.5 Rads Per Sec 



Sea State 
Kl = 

Tl= 

T2 = 

T3 = 

T4 = 

COST J= 



2.3292980 
0.5873335 
21.2806244 
20.8198700 
5.3823652 
0. 1791609E-02 



2.3723307 
0 . 5848479 
8.7543087 
7 . 1684380 
7.0838490 
0 . 7 157017E-02 



2.3858051 
0.6092799 
8.3835888 
6.6424551 
7.5094414 
. 2188155E-01 
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TABLE 11 



Controller C For 90° Encounter Wave Angle 



Encounter Frequency 0.2 Rads Per Sec 



Sea State 
Kl= 

Tl = 

T2 = 

T3 = 

T4 = 

COST J= 



Sea State: 
Kl= 

Tl = 

T2 = 

T3 = 

T4= 

COST J= 



1.3517475 
35.3476257 
19.6020966 
15.6134491 
19 . 6564178 
0 . 2690107E+01 



1.3857212 
36.2312164 
21.4780426 
15.0542984 
21.3916779 
0. 9237458E+01 



1.4623213 
42. 9548187 
28.4767914 
14.1265106 
28.5155182 
0 . 2028659E+02 



Encounter Frequency 0.4 Rads Per Sec 



4 

1.0069885 
49.7900391 
36.0238495 
15.7164612 
35.5061340 
0 . 8278778E + 01 



1.1330795 
3.6114454 
18.1521301 
31.1608582 
17.9342651 
0 . 1712465E+02 



1 

0 

37 
39 

38 



1926146 

7214398 

0983887 

9669037 

7577972 



0 . 2298 665E + 02 



Encounter Frequency 0.6 Rads Per Sec 



4 

1.8570919 
1.0337286 
11.9529600 
12.3197098 
11.7377634 
COST J= 0 . 9655529E-01 



Sea state: 
Kl= 

Tl = 

T2 = 

T3 = 

T4 = 



1.9276266 
0.9287271 
15.0574036 
13.4190674 
15.0851593 
0 . 2932869E+00 



1.8934374 
0.4949190 
30.3212738 
27.2785339 
26.5594635 
0 . 5224001E + 00 



Encounter Frequency 0.75 Rads Per Sec 



Sea state: 
Kl= 

Tl = 

T2 = 

T3 = 

T4 = 



4 

1.9137421 

0.4686716 

26.6723328 

20.3362122 

30.0011749 



1 

0 

32 

27 

32 



8589468 

3916095 

9264069 

4372711 

4602356 



COST J= 0 . 4419509E-01 0 . 1722615E+00 



1.7544203 
0. 1022463 
42.4169006 
29 . 9568634 
52.5700073 
0 . 4958326E+00 



Encounter Frequency 1.5 Rads Per Sec 



Sea State 
Kl= 

Tl= 

T2 = 

T3 = 

T4 = 

COST J= 



4 

2.3951197 
0.6277218 
4.0210686 
6.0798759 
69.4311066 
0. 1650761E-01 



2.4414005 
0.6339726 
6.7097816 
4.9043417 
6.7554045 
0. 6472868E-01 



2.4639597 
0.6183524 
7.3951139 
5 . 1939421 
7 . 3074417 
0. 1906425E+00 
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Figure 4.2 Yaw Vs Time, Sea State 4. 
Encounter Frequency 0.2 rads per sec , Encounter Angle 30 
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Figure 4.3 Yaw Vs Time, Sea State 4. 
Encounter Frequency 0.4 rads per sec , Encounter Angle 30 
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Figure 4.4 Yaw Vs Time, Sea State 4. 

Encounter Frequency 0.6 rads per sec, Encounter Angle 30 
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Figure 4.6 Yaw Vs Time, Sea State 4. 

Encounter Frequency 1.5 rads per sec, Encounter Angle 30 
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Figure 4.7 Yaw Vs Time, Sea State 4. 

Encounter Frequency 0.2 rads per sec, Encounter Angle 60 
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Figure 4.8 Yaw Vs Time, Sea State 4. 

Encounter Frequency 0.2 rads per sec, Encounter Angle 90 




38 



Figure 4.9 Rudder Vs Time, Sea State 4. 
Encounter Frequency 0.2 rads per sec, Encounter Angle 30 
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Figure 4.10 Rudder Vs Time, Sea State 4. 
Encounter Frequency 0.4 rads per sec, Encounter Angle 30 
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Figure 4.11 Rudder Vs Time, Sea State 4. 
Encounter Frequency 0.6 rads per sec. Encounter Angle 30 
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Figure 4.12 Rudder Vs Time, Sea State 4. 
Encounter Frequency 0.75 rads per sec, Encounter Angle 30 
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Figure 4.13 Rudder Vs Time, Sea State 4. 
Encounter Frequency 1.5 rads per sec, Encounter Angle 30 
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Figure 4.14 Rudder Vs Time, Sea State 4. 
Encounter Frequency 0.2 rads per sec. Encounter Angle 60 
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Figure 4.15 Rudder Vs Time, Sea State 4. 
Encounter Frequency 0.2 rads per sec, Encounter Angle 90 
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Figure 4.16 Cost Vs Tl, Sea State 7. 

Encounter Frequency 0.2 rads per sec, Encounter Angle 30 
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Figure 4.17 Cost Vs T2 , Sea State 7. 

Encounter Frequency 0.2 rads per sec, Encounter Angle 30 
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Figure 4.18 Cost Vs T3 , Sea State 7. 

Encounter Frequency 0.2 rads per sec. Encounter Angle 30 
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Figure 4.19 Cost Vs T4 , Sea State 7. 

Encounter Frequency 0.2 rads per sec, Encounter Angle 30 
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Figure 4.20 Cost Vs Kl, Sea State 7. 

Encounter Frequency 0.2 rads per sec, Encounter Angle 30 



V. IRREGULAR SEAS -CONTROLLER DESIGN 



A sea state generator program which genarates added 
mass, added inertia values and in addition calculates the 
forces and moments, was coupled to the fortran program shown 
in Appendix D, so that the function minimization subroutine 
(BOXPLX) could be used in the presence of the irregular 
sea. The forces and moments were stored in a look up table 
which was coupled to the equation of motion. 

The optimal gains obtained for 0.75 rad/sec encounter 
frequency in the regular sea study were used as the initial 
guess in order to evaluate the optimal controller parameters 
involving the irregular sea. 

It is known that sea is never regular but actually is a 
random phenomenon where waves are continually changing in 
height, length and breadth. 

Since the sea state during this study is represented by 
irregular waves then the waves impinding on the ship hull 
would contain the total energy density spectrum. This 
dencity specrum is composed of many frequencies and there- 
fore the response of the ship would be to an average value 
of added mass and added inertia. 

For this study where the ship's speed is 23 knots, the 
controller has the form of equation B.l, we used values for 
added mass and added inertia corresponding to 0.75 rad/ sec 
because close to that frequency the energy density is 
maximum. 

The optimal controller parameters found are shown in 
Table 12, and typical system's response is shown in Figures 
5. 1,5. 2, 5. 3, 5. 4, 5. 5 

These plots were obtained using the program of Appendix 
E. 
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A careful analysis of the extracted results leads to 
conclude : 

• The maximum deviation of controller parameters values 
occured at 30° encounter angle. 

• For all encounter angles the maximum cost occured for 
sea state 7. 

• For specific encounter angles the higher the sea state 
the higher the cost. 
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TABLE 12 

Controller C for Different Encounter Wave Angle 

0° Encounter Wave Angle 
Encounter Frequency 0.75 Rads Per Sec 



Sea State: 
Kl = 

Tl= 

T2 = 

T3 = 

T4= 

Cost J= 



1.7460003 
35 . 2969971 
22.2480011 
13 . 6510000 
22.1170044 
0 . 9438775E- 34 



6 

1.7460003 
35.2969971 
22.2480011 
13.6510000 
22.1170044 
0 . 3782326E-34 



0 



1.7460003 
35.2969971 
22.2480011 
13 . 6510000 
22.1170044 
4605412E- 33 



30° Encounter Wave Angle 
Encounter Frequency 0.75 Rads Per Sec 



Sea State: 


4 


6 


7 


Kl = 


2.4596768 


1.5688477 


2.4804857 


Tl= 


88.2797241 


47.3260040 


56 . 3383179 


T2 = 


50.5678864 


35 . 6789203 


51.4950714 


T3 = 


5.2703905 


21.5429993 


5 . 7071409 


T4= 


95 .3189392 


25.0237122 


91.6153102 


Cost J= 0. 


8905333E-01 


0 . 4461777E-01 


0. 1304857E+00 



60° Encounter Wave Angle 
Encounter Frequency 0.75 Rads Per Sec 



Sea State: 


4 


6 


7 


Kl= 


0.3419376 


0.1919854 


0.2646747 


Tl= 


99.6322021 


16.3320312 


15.4308896 


T2 = 


35.0456085 


38.0289001 


35.0908203 


T3 = 


38 . 9945677 


34.9145813 


38.4357300 


T4= 


25 . 1149750 


60.8223572 


62.9836578 


Cost J= 0. 


1788034E-03 


0 . 3666982E-02 


0 . 1193313E-01 



90° Encounter Wave Angle 
Encounter Frequency 0.75 Rads Per Sec 



Sea State 
Kl= 

Tl= 

T2 = 

T3 = 

T4= 

Cost J= 



4 

0.5902819 
87.4059400 
34.5666444 
39.6794654 
25.1668799 
0. 7864684E-01 



2.4912267 
55.7961656 
35.6789999 
30.2094433 
25.0789005 
0. 1778146E-01 



2 

32 

33 
4 

25 



,4547853 

,7171021 

8865433 

9768888 

3325553 



0 . 2522047E + 00 
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Figure 5.5 Rudder vs Time, Sea State 6. 
Encounter Frequency 0.75 Rads Per Sec , Encounter Angle 30 



VI. ADAPTIVE CONTROL 



A. STRUCTURE 

When the ship is moving in a seaway, the controller 
parameters are changing due to alterations in the sea state 
and encounter angle. In addition we know that using fixed 
parameters for the controller over the entire spectrum, it is 
somewhat difficult to have an appropriate response of the 
system. The adjustment of the controller parameters during 
operation in a seaway, can be achieved by means of an adap- 
tive control. 

The adaptive controller can be built with digital 
circuits and analog components as well. Analog system hard- 
ware has to be designed for each specific requirement, and 
any new requirement involves changes to components. This is 
a time consuming task. However , for simple control require- 
ments anolog systems still have a possible economic advan- 
tage over digital systems [Ref. 11]. 

On the other hand, digital systems are immediately more 
attractive when control systems are required to carry out 
more and more complex tasks. The advent of microprocessors 
and associated components has enabled low cost microcom- 
puters to be built. These no longer require special environ- 
ments and are fully compatible with shipboard use. The 
advantages of microcomputer controls over other types are 
that they are extremely reliable particularly as relatively 
fewer components are required and hence they are smaller. 
Their capability is greater than comparable analog systems 
due to their ability to carry out more complex calculations. 
A major advantage is their flexibility while using standard 
hardware, being able to be reconfigured for changes in 
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system requirement without the need to alter the hard- 
ware. This flexibility is achieved by the programmed software 
which is stored in the memory of the computer. 

An adaptive control scheme is indicated in Figure 6.1 
and in Figure 6.2 we show analytically the components of 
what we call the 'Decision Generator'. 

B. SEARCH ALGORITHM 

The parameters corresponding to various environmental 
conditions (sailing modes) are sorted in a memory. When the 



value of the performance index changes from its theoretical 
value. This change is detected by the threshold detector. 

If the threshold is too small the adaptive structure 
will try new conditions when it is not necessary and in the 
opposite case the adaptive structure might not adapt to new 
sailing modes. 

Using the term threshold we want to make sure that only 
significant changes in the environmental conditions will be 
considered, otherwise the system will look for a new model 
too often. 

When the threshold detects that a new set of parameters 
should be used the search algorithm will look for conditions 
close to the previous one because the external conditions 
are not expected to vary rapidly. Then a new optimal value 
of J is selected and must be matched with the performance 
index computed for the real psi and delta. 

Matching operation must be done after some delay 
ensuring that the conditions had an effect on the real cost 
function. 

These new parameter values will be fed into the 
controller and function minimization subroutine as well. 



sailing modes change, 




and psi(^p) vary so that the 
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Figure 6.1 Adaptive Control Scheme 
C. ADVANTAGES 

• Using this adaptive scheme we don't use sensors which 
may be unrealistic and the use of predetermined filters 
which are so expensive can be avoided. 
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Figure 6.2 Decision Generator 

• Provides a good choice for the function minimization 
subroutine to start it's iteration algorithm and there- 
fore we get the optimum values as quickly as possible. 

D. DISADVANTAGES 

• The search algorithm must be carefully determined. 

• Since we don't include sensor's measurements there is 
no indications on how to perform the search. This 
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problem can be eliminated when the NAVSTAR/ GLOBAL 
POSITION SYSTEM (GPS) will be able to provide precision 
navigation data. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

• A properly designed controller can minimize the rudder 
activity providing desired overseas heading as well, and 
therefore does result in substantial improvement in 
propulsion efficiency. 

• Actual savings in fuel cannot yet be determined since 
there is no information available from the conventional 
autopilot and therefore there is no possibility for 
comparison . 

• It is believed that the performance index used in this 
research is a fairly adequate function. Doubts arise 
from the weighting factor which is included and this 
because lambda is based on assumptions and it's accu- 
racy is not certain. 

• An adaptive controller which minimizes propulsion 
losses due to steering is needed when environmental 
conditions and ship characteristics change. 

• Studying all the investigated sailing modes, it turns 
out that the cost surface is flat. As a consequence 
determination of the controller parameters does not 
require high accuracy. 

• The study shows that the use of the third order ship 
Nomoto model is a reasonable choice instead of using 
the ship's equation of motion, which involves both the 
sway and yaw equations. 

• It can be assumed from the work done in this research 
that some of the findings could have applicability to a 
number of ship types . However it is not possible to make 
inference, of a general nature concerning all ships 
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from the results of the SL-7 which represents a partic- 
ular type of a particular class of ships. 

B. RECOMMENDATIONS FOR FUTURE STUDY 

• Additional work has to be done for obtaining optimum 
controller parameters under an expanded range of oper- 
ating conditions. The more optimum controller param- 
eter values available, the better the determination of 
the search algorithm recommended in the adaptive 
control scheme. 

• A study including the surge equation in our ship model 
is recommended since in reality added resistance due to 
steering reduces ship's speed and so far we assumed 
constant speed. In addition with the use of the surge 
equation we can determine actual energy losses. 

• Investigation of more advanced methods of adaptive 
control based on "on-line" determination of plant 
characteristics . 

• So far we were interested in minimizing rudder and 

yawing activity to reduce propulsion losses using a 
particular performance criterion in open 

seas . Considering some other types of control such as 
automatic control for replenishment requires an inves- 
tigation as to the suitability of the cost function, 
and a comparison with other potential cost functions. 
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APPENDIX A 



DETERMINATION OF NOMOTO THIRD ORDER MODEL-BOXPLX 

//PROGRA JOB (????, 0356) , 'RESEARCH' ,CLASS=G 

//“MAIN 0RG=NPGVM1. ????P 

// EXEC FORTXCG , PARM . FORT= ' OPT ( 2 ) ' , IMSL=DP , REGION= 1024K 

/ /FORT . SYSIN DD * 

C THIS PROGRAM WILL OBTAIN THE CONTROLLER OPTIMAL GAINS. 

C IT IS REFERENCED IN CHAPTER 5. 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 

C OBTAINED CHANGE XS(*) TO X(*) AND DELETE XU(*),AND XL(*). 
DIMENSION X S ( 5 ) , XU ( 5 ) , XL ( 5 ) 

XS(1)=1. 7466650 
XS(2) = 35. 2978973 
XS(3) = 22. 2485657 
XS(4)=13. 6511078 
XS(5)=22. 1172638 

C XS(I) IS THE STARTING GUESS 

C XL(I) IS THE LOWER LIMIT FOR THE I ' TH VARIABLE 

C XU (I) IS THE UPPER LIMIT FOR THE I ' TH VARIABLE 
XL ( 1 ) = 0 . 1 
XU ( 1 ) = 4 . 0 
XL ( 2 ) = 0 . 1 
XU(2)=80.0 
XL ( 3 ) = 1 . 0 
XU(3)=50.0 
XL (4) =1.0 
XU(4)=30.0 
XL ( 5 ) = 1 . 0 
XU(5)=50.0 

C A DESCRIPTION OF THE FOLLOWING PARAMETERS 

C IS DISCUSSED IN BOXPLX 
R= 9 . / 13 . 
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NTA= 1000 
NPR= 100 
NAV=0 
NV=5 
IP=0 

C THE FOLLOWING STATEMENT MUST BE CHANGED TO 
C CALL PLANT (X) 

C IF ONLY SIMULATION IS WANTED 

CALL BOXPLX (NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER) 
WRITE (6,25) 

25 FORMAT (IX,' OPTIMAL GAINS',/) 

DO 30 1=1,5 

30 WRITE (6, 40)1, XS (I) 

40 FORMAT ( IX , ' X( ' , 12 , ' ) = ' , F14 . 7 ) 

STOP 

END 

SUBROUTINE PLANT (XX) 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
COMMON TDIFF 
REAL* 8' L,L2,L3,L4,L5,L6 

REAL* 8 X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL*8 TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL* 8 YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
REAL* 8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL*8 RHO , IZ , FX , FY ,MZ ,XP ,MASS , DELT ,MZI , RXI , WA , WE 
REAL* 8 DYAWE , YAWE , YAWC , I SE , I SR , LAMDA , D , RYR , RYI , MZR 
REAL*8 K1 , T1 ,T2 ,D ,X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR 
REAL*8 T3 ,T4 ,X3 ,DX3 ,X4 
DIMENSION XX (5 ) 

C 

C CLOSE LOOP ANALYSIS WITH FILTER 
C 

C INITIAL CONDITIONS FOR INTEGRATION 
C SIMULATION END TIME IN SECONDS 
ETIME= 600.0 



66 



TIME=0 . 0 
ICOUNT= 1 

C INITIALIZE THE COST FUNCTION 
ISE=0 . 0 
ISR=0 . 0 
TDIFF=0 . 0 
LAMDA=8. 128 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
K1=XX( 1) 

T1=XX(2 ) 

T2=XX( 3 ) 

T3=XX(4) 

T4=XX(5) 

C WRITE(6 , 1010) K1,T1,T2 

C1010 FORMAT ( IX , 'K1 =',F15.7,'T1 =',F15.7, ? T2 =',F15.7) 

C X , XDOT , Y , YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0 . 0 
XDOT=0 .0 
YDOT=0 . 0 

C U , UDOT , V , VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT=0 .0 
VDOT=0 . 0 
YAW=0 . 0 
R=0 .0 
RDOT=0 . 0 

C ORDERED SPEED IN FEET/ SEC 

C 38.82 FT/ SEC=23 KNOTS 
UC= 38 . 82 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 

C D = RUDDER ANGLE 
D=0 . 0 
L=880 . 5 
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L2=L**2 

L3=L*L*L 

L4=L- V L3 

L5=L*L4 

L6=L*L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
C MOMENTS IN Z 
FX=0 . 

FY=0 . 

MZ = 0 . 

RXR= - 0 . 15744D+05 
RXI=-0 . 19950D+06 
RYR=0 . 52365D+04 
RYI=0 . 18699D+06 
MZR=-0 .29870D+08 
MZI=-0 . 35751D+07 
C RXR=-0 .50230D+04 

C RXI=0 . 12712D+05 

C RYR=0 . 35290D+04 

C RYI=-0 .31909D+05 

C MZR=0 .38826D+07 

C MZI=-0 .643I3D+07 

C RXR=0 . 28540D+04 

C RXI=-0 .99574D+04 

C RYR= -0 . 85441D+04 

C RYI=0 .39595D+05 

C MZR=-0 . 13014D+08 

C MZI=0 . 11348D+08 

C RXR=-0 . 75642D+04 

C RXI=0 .83497D+04 

C RYR=0 . 23379 D+ 05 

C RYI=-0 .8150 2D +05 

C MZR=0 .28622D+07 

C MZI=-0 . 19388D+08 
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c 



c 



c 



c 



c 



c 



RXR= -0.37916D+04 
RXI = 0 . 1638 1D+ 04 
RYR= -0.76647D+05 
RYI= -0. 37685D+05 
MZR= -0.83915D+07 
MZI=-0. 53176 D+ 07 



RX=DSQRT (RXR**2+RXI**2 ) 

RY= DSQRT ( RYR**2 + RYI**2 ) 

RZ=DSQRT (MZR**2+MZI**2 ) 

TX= DAT AN ( RXI / RXR ) 

T Y = DAT AN ( RYI / RYR ) 

TZ=DATAN(MZI/MZR) 

C SIGNIFICANT WAVE HEIGHT; SEA STATE 1-5,2-10,3-15 

C 4-17.5,5-22.5,6-27,7-35,8-42,9-60 
WA= 17.5 

C ENCOUNTER FREQUENCY . 1 , . 2 , . 3 , . 4 , . 5 , . 6 , . 7 5 , 1 . 0 , 1 . 5 , 2 . 5 



C HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE AS 
c PARAMETERS 

RHO=l. 9876 

MASS= ( .0044)*( .5*RHO*L3) 

IZ= (0 . 00028 )*( . 5"RHO“L5 ) 

YAWE=0 . 0 
X2=0 . 0 
DX2=0 . 0 
X3 = 0 . 0 
DX3=0 .0 
200 CONTINUE 

S=DSQRT (u**2+V**2 ) 

C INPUT YAW COMMAND 
YAWC=0 . 0 

IF (TIME. GE. 0.0) YAWC=0.0 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW ORDERED) 
C ( COMPENSATOR FILTER ) 

YAWE=YAW - YAWC 



WE = 0 . 2 
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DX2 = (YAWE-X2 ) /T2 
X4=K1*(T1*DX2+X2) 

DX3= (X4-X3 ) /T4 
D= (T3“DX3+X4) 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

C XU DOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , 

C ENCOUNTER FREQUENCY 

XUDOT = ( - . 0 0 0 1 ) * ( . 5 *RHO*L3 ) 

XU= ( - 0 . 025 3 ) * ( . 5 *RHO*L2*S ) 

XUU= ( - 0 . 0 0 0 3 ) * ( . 5*RHO*L2 ) 

XVR= ( 0 . 0 0 3 9 ) * ( . 5 "RHO"L 3 ) 

XVV= ( - . 0012 )*( . 5*RHO*L2 ) 

XDD= ( -0 . 0005 )*( . 5 *RHO *L2 * S * * 2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

C YV= ( - 0 . 00 7 5 8 ) * ( . 5 *RHO*L2*S ) 

YR= ( 0 . 00 2 3 ) * ( . 5 *RHO*L3 *S ) 

YD= (0 . 00145 )*( . 5*RHO*L2*S**2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= ( - 0 . 0 0 8 ) * ( . 5 *RH0*L4 / S ) 

YVVV= ( - 0 . 0 3 ) * ( . 5 *RHO*L2 / S ) 

YRRR= ( 0 . 0 0 3 ) * ( . 5 *RHO*L5 / S ) 

YDDD= ( -0 . 0005 )*( . 5*RHO*L2*S**2 ) 

C YUDOT IS THE ADDED MASS TERM WHICH MUST BE 
C CHANGED FOR DIFFERENT ENCOUNTER ANGLE , SPEED , 

C ENCOUNTER FREQUENCY 
C YVDOT= ( - 0 . 003 9 ) * ( . 5*RHO*L3 ) 

C SPEED=23 KNOTS, ENCOUNTER ANGLE = 30 , ENCOUNTER 
C FREQUENCY =0.2 

YVDOT= -0.30908D+07 
YV= -0.8127 1D+04 
C YVDOT= -0.36185D+07 

C YV= -0.24757D+06 

C YVDOT = -0.32890D+07 

C YV=-0.11775D+07 
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C YVDOT= - 0.23038D+07 

C YV=-0. 18267D+07 

C YVDOT= - 0.59800D+06 

C YV=-0. 13260D+07 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV= ( - 0 . 00213 ) * ( . 5*RHO*L3*S ) 

C NR= ( - 0 . 00 105 ) * ( . 5*RHO*L4*S ) 

ND= ( - 0 . 0007 )* ( . 5*RHO*L3*S**2 ) 

NVVR= ( - 0 . 015 )*( . 5*RHO*L4 / S ) 

NVRR= ( - 0 . 008 ) * ( . 5 *RHO*L5 / S ) 

NVVV= ( 0 . 0 1 ) * ( . 5*RHO*L3 / S ) 

NRRR= ( - 0 . 006 ) * ( . 5*RHO*L6 / S ) 

NDDD= (0 . 0001 )*( . 5*RHO*L3*S**2 ) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE 
C CHANGED FOR DIFFERENT ENCOUNTER ANGLE , SPEED , 

C ENCOUNTER FREQUENCY 
C 

C NRDOT = ( -0 . 00027 )*( . 5*RHO*L5 ) 

C SPEED=23 KNOTS, ENCOUNTER ANGLE = 90 , ENCOUNTER 
C FREQUENCY =0.2 

NRDOT=-0 . 26251D+ 12 
NR= -0.53637D+09 
C NRDOT= -0 . 20125D+ 12 

C NR= -0 . 94970D+ 10 

C NRDOT= -0. 18671D+12 

C NR= -0.46860D+11 

C NRD0T=-0 . 14518D+12 

C NR= -0 .87538D+11 

C NRDOT= -0. 37261D+11 

C NR= -0 .69856D+11 

C REGULAR WAVE SEA STATE 

FX=WA*RX*DCOS (WE*TIME+TX) 

F Y = WA*RY*DCO S ( WE*T IME + T Y ) 

MZ = WA*RZ*DCO S ( WE*T IME + T Z ) 

C U ACTUAL SPEED 
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C UC COMMANDED SPEED 
C XP = PROPELLER THRUST 
XP=-XUU*UC**2 
C EQUATIONS OF MOTION 

C UDOT= ( (XVR + MASS )*V*R + XUU*U**2 + XVV*V**2 

C 1 + XDD*D*D + FX + XP ) / (MASS-XUDOT ) 

VDOT= (YV*V + ( YR- MAS S*U ) *R + YD-'D + YVVR*V**2*R 

1 + YVRR*V*R**2 +YRRR*R**3 

2 + YDDD*D**3 + FY ) / (MASS- YVDOT) 

RDOT= ( NV*V + NR*R + ND*D + NVVR*V**2*R 

1 + NVRR*V*R**2 +NVW*V**3 

2 + NRRR*R**3 + NDDD*D**3 + MZ )/(IZ-NRDOT) 

C WHEN TO PRINTOUT 

IF (ICOUNT.EQ.il) GO TO 50 
GO TO 300 

C CONVERT RADIANS TO DEGREES 
5 0 YAWDEG = YAW* 57.296 
RDEG = R*5 7.296 
RDDEG=RDOT*57 .296 
DDEG=D*57 .296 
YAWC=YAWC*57 .296 

C WRITE (6,100) TIME , XP , X , XDOT , Y , YDOT 

C 1 ,UC,U,UDOT,V,VDOT,YAWC, YAWDEG, RDEG, RDDEG,DDEG 
100 FORMAT ( IX , ’ TIME= ’ , F8 . 3 , ’ SEC XP= ? ,FI0.2,’ LBF 
1 X=’ ,F8.2, 

1’ FT XDOT= ’ , F8 . 4 , ’ FT/ SEC Y=’,F8.2,’ FT YDOT 
1 , F8 . 4 , ’ FT/SEC 

1 ' , / , 2X , ' UC= ' , F8 . 4 , ' FT/SEC U=',F8.4,' FT/SEC 
1 UDOT= ' , F10 . 6 , 

1 ’ FT / SEC**2 V= ' , F8 . 4 , ' FT/ SEC VDOT=’,F10.6 

1 ' FT/SEC**2 ' 

1 , / ,2X, ’YAWC= ’ ,F8 .4, ' DEG YAW= ' 

1 , F15 . 7 , ' DEG YAW RATE= ' , F15 . 7 , ' DEG/ SEC 

1 YAW ACCEL= ' , F15 . 7 , ' 

1 DEG / SEC**2' ,/ ,2X, 'RUDDER =’,F15.7,' DEG ',/) 
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ICOUNT= 1 

C TEST IF WANT TO STOP 
300 IF (TIME . GE . ETIME ) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELT= 1 . 0 
C INTEGRATION 

U=U+UDOT*DELT 
V=V+VDOT"DELT 
R=R+RDOT*DELT 
YAW = YAW + R " DELT 
X2=X2+DX2*DELT 
X3 = X3 + DX3 -'DELT 

C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
C XDOT=U"DCOS (YAW) -V*DSIN (YAW) 

C YDOT=U "DSIN ( YAW ) + V*DCOS ( YAW ) 

C X=X+XDOT*DELT 

C Y=Y+YDOT"DELT 

TIME= TIME + DELT 
ICOUNT=ICOUNT+ 1 
ISE = ISE + L AMD A -Y AWE * * 2 
ISR= ISR + D**2 
GO TO 200 

C J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 

WRITE (6, 500) ISE , I SR , TDIFF ,Kl,Tl,T2,T3 ,T4 
500 FORMAT (’ ' , IX , ’ TOTAL= ’ , F15 . 7 , 2X , 

1 'Kl=' ,F15.7,2X, 'Tl=' ,FI5.7,2X, 'T2=' ,F15.7,2X, 

1 'T3=' ,F15.7,2X, ’ T4= ' , FI5 . 7) 

RETURN 

END 

SUBROUTINE BOXPLX (NV ,NAV ,NPR, NTZ ,RZ ,XS , IP , BU , BL , 
1YMN , IER) 

C 

DIMENSION V(50, 50) , FUN(50), SUM(25), CEN(25), 
1XS(NV) ,BU(NV) , BL (NV ) 
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c 



KV = 5 
EP = l.E-6 
NTA = 2000 

IF (NTZ.GT.O) NTA = NTZ 
R = RZ 

IF (R.LE.O. .OR.R.GE. 1. ) R=l./3. 

NVT = NV+NAV 
C 

C TOTAL VARS, EXPLICIT PLUS IMPLICIT 

NT = 0 

C CURRENT TRIAL NO. 

NPT = 0 

C CURRENT NO. OF PERMISSIBLE TRIALS 

NTFS = 0 

C CURRENT NO. OF TIMES F HAS BEEN ALMOST UNCHANGED 

C 

C CHECK FEASIBILITY OF START POINT 

C 

DO 4 1=1, NV 
VT = XS(I) 

IF (BL(I) .LE.VT) GO TO 1 
II = -I 
VT = BL(I) 

GO TO 2 

1 IF (BU(I) .GE.VT) GO TO 3 
II = I 

VT = BU ( I ) 

2 IF (NPR.GT.O) WRITE (6,49) II 

3 V(I,1) = VT 
CEN(I) = VT 

IF (IP.EQ.l) GO TO 4 

BL(I) = BL(I)+AMAX1(EP,EP*ABS(BL(I) ) ) 

BU ( I ) = BU ( I ) -AMAX1(EP , EP'vABS (BU (I ) ) ) 

4 SUM ( I ) = VT 
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c 

c 

NCE = 1 

C NUMBER OF CONSTRAINT EVALUATIONS 

1=1 

IF (KE(V (1,1)). EQ . 0 ) GO TO 5 
IF (NPR.LE.O) GO TO 12 
WRITE (6,50) 

GO TO 12 

5 NFE = 1 
C 

C NUMBER OF VERTICES (K) = 2 TIMES NO. OF VARIABLES. 

K = 2*NV 
C 

C NUMBER OF DISPLACEMENTS ALLOWED. 

NLIM = 5*NV+10 
C 

C NUMBER OF CONSECUTIVE TRIALS WITH UNCHANGED 
c FE TO TERMINATE. 

NCT = NLIM+NV 
ALPHA = 1.3 
FK = K 
FKM = FK-1. 

BETA = ALPHA + 1 . 

C 

C INSURE SEED OF RANDOM NUMBER GENERATOR IS ODD. 

IQR = R*1 . E7 

IF (MOD (IQR, 2) .EQ.O) IQR=IQR+101 
C 

C SET UP INITIAL VERTICES 

FUN ( 1 ) = FE(V ( 1 , 1 ) ) 

YMN = FUN ( 1 ) 

6 FI = 1. 

FUNOLD = FUN ( 1 ) 

C 
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DO 15 1=2, K 
FI = FI+ 1 . 

LIMT = 0 

7 LIMT = LIMT+1 
C 

C END CALCULATION IF FEASIBLE CENTROID CANNOT BE FOUND. 

IF (LIMT.GE.NLIM) GO TO 11 
C 

DO 8 J = 1 , NV 
C 

C RANDOM NUMBER GENERATOR (RANDU) 

IQR = IQR*65539 

IF (IQR. LT . 0 ) IQR = IQR+2147483647+ 1 
RQX = IQR 

RQX = RQX* . 4656613E- 9 

V(J,I) = BL ( J ) + RQX* (BU(J)-BL(J)) 

IF (IP.EQ.l) V( J , I ) =AINT (V ( J , I ) + .5) 

8 CONTINUE 
C 

DO 10 L= 1 , NLIM 
NCE = NCE+1 

IF (KE (V ( 1 , I ) ) . EQ . 0 ) GO TO 13 
C 

DO 9 J= 1 , NV 

VT = . 5*(V( J , I ) +CEN ( J) ) 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 

V(J,I) = VT 

9 CONTINUE 
C 

10 CONTINUE 
C 

11 IF (NPR.LE.O) GO TO 12 
WRITE (6,51) I 

CALL BOUT (NT , NPT , NFE , NCE , NV , NVT , V , I , FUN , CEN , I ) 

12 IER = -1 
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GO TO 48 



C 

13 DO 14 J = 1 , NV 

SUM(J) = SUM( J ) +V ( J , I ) 

14 CEN(J) = SUM(J)/FI 
C 

C TRY TO ASSURE FEASIBLE CENTROID FOR STARTING. 

NCE = NCE+1 

IF (KE(CEN) .EQ.O) GO TO 60 
SUM ( J ) = SUM(J) -V(J,I) 

GO TO 7 

60 NFE = NFE+1 

FUN ( I ) = FE (V( 1 , I ) ) 

15 CONTINUE 
C 

C END OF LOOP SETTING OF INITIAL COMPLEX. 

IF (NPR.LE.O) GO TO 17 

CALL BOUT (NT , NPT , NFE , NCE , NV , NVT , V ,K , FUN , CEN , 0 ) 
C 

C FIND THE' WORST VERTEX, THE ’J’TH. 

J = 1 
C 

DO 16 1=2, K 

IF ( FUN ( J ) . GE . FUN ( I ) ) GO TO 16 
J = I 

16 CONTINUE 
C 

C BASIC LOOP. ELIMINATE EACH WORST VERTEX 
C IN TURN. IT MUST BECOME NO LONGER WORST, NOT 
C MERELY IMPROVED. FIND NEXT-TO-WORST VERTEX, 

C THE 'JN’TH ONE. 

17 JN = 1 

IF (J.EQ. 1) JN = 2 
C 

DO 18 1=1, K 
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IF (I.EQ.J) GO TO 18 
IF (FUN(JN) .GE.FUN(I) ) GO TO 18 
JN = I 

18 CONTINUE 
C 

C LIMT = NUMBER OF MOVES DURING THIS TRIAL TOWARD 

C THE CENTRIOD DUE TO FUNCTION VALUE. 

LIMT = 1 
C 

C COMPUTE CENTROID AND OVER REFLECT WORST VERTEX. 
C 

DO 19 1=1, NV 
VT = V(I,J) 

SUM ( I ) = SUM( I ) - VT 

CEN(I) = SUM ( I ) / FKM 

VT = BETA*CEN(I)-ALPHA*VT 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 

C 

C INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED. 

19 V(I,J) = AMAX1 (AMIN1 (VT , BU ( I ) ) , BL( I ) ) 

C 

NT = NT+ 1 
C 

C CHECK FOR IMPLICIT CONSTRAINT VIOLATION. 

C 

20 DO 25 N= 1 , NLIM 
NCE = NCE+1 

IF (KE (V ( 1 , J ) ) . EQ . 0 ) GO TO 26 
C 

C EVERY 'KV'TH TIME, OVER- REFLECT THE OFFENDING 
C VERTEX THROUGH THE BEST VERTEX. 

IF ( MOD ( N , KV ) . NE . 0 ) GO TO 22 
CALL FBV ( K , FUN , M ) 

C 

DO 21 1=1, NV 
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VT = BETA'' V ( I , M ) - ALPHA "V ( I , J ) 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 

21 V(I,J) = AMAX1 ( AMINl ( VT , BU (I ) ) , BL ( I ) ) 

C 

GO TO 24 
C 

C CONSTRAINT VIOLATION: MOVE NEW POINT TOWARD CENTROID 

C 

22 DO 23 1=1, NV 

VT = . 5 * ( CEN (I)+V(I,J)) 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 

V(I,J) = VT 

23 CONTINUE 
C 

24 NT = NT+ 1 

25 CONTINUE 
C 

IER = I 
C 

C CANNOT GET FEASIBLE VERTEX BY MOVING TOWARD CENTROID, 
C OR BY OVER-REFLECTING THRU THE BEST VERTEX. 

IF (NPR.LE.O) GO TO 42 
WRITE (6,52) NT , J 

CALL BOUT (NT , NPT ,NFE , NCE , NV ,NVT , V ,K , FUN , CEN , J ) 

GO TO 42 
C 

C FEASIBLE VERTEX FOUND , EVALUATE THE OBJECTIVE FUNCTION 

26 NFE = NFE+1 
FUNTRY = FE (V ( 1 , J ) ) 

C 

C TEST TO SEE IF FUNCTION VALUE HAS NOT CHANGED. 

AFO = ABS ( FUNTRY -FUNOLD) 

AMX = AMAX1 (ABS (EP '"FUNOLD ) ,EP) 

C 

C ACTIVATE THE FOLLOWING TWO STATEMENTS 
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C FOR DIAGNOSTICS PURPOSES ONLY. 

C WRITE (6,99) J , AFO , AMX , FUNTRY , FUNOLD , FUN ( J ) , 

C 1FUN(JN) ,NTFS,N 
C 99 FORMAT ( IX , 13 , 6E15 . 7 , 215 ) 

IF (AFO. GT. AMX) GO TO 27 
NTFS = NTFS+ 1 
IF (NTFS.LT.NCT) GO TO 28 
IER = 0 

IF (NPR.LE.O) GO TO 42 
WRITE (6,53) K 

CALL BOUT (NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , 0 ) 

GO TO 42 

27 NTFS = 0 
C 

C IS THE NEW VERTEX NO LONGER WORST? 

28 IF (FUNTRY.LT. FUN ( JN ) ) GO TO 34 
C 

C TRIAL VERTEX IS STILL WORST; ADJUST TOWARD CENTROID. 

C EVERY 'KV'TH TIME, OVER- REFLECT THE OFFENDING 

C VERTEX THROUGH THE BEST VERTEX. 

LIMT = LIMT+1 

IF (MOD (LIMT ,KV ) . NE . 0 ) GO TO 30 
CALL FBV (K , FUN ,M) 

C 

DO 29 1=1, NV 

VT = BETA - v V (I,M)- ALPHA - v V (I, J) 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 

29 V(I,J) = AMAX1 ( AMIN1 (VT , BU ( I ) ) , BL ( I ) ) 

C 

GO TO 32 
C 

30 DO 31 1=1, NV 

VT = . 5* ( CEN (I)+V(I,J)) 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 

V(I,J) = VT 
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31 CONTINUE 



C 

32 IF (LIMT.LT.NLIM) GO TO 33 
C 

C CANNOT MAKE THE ' J ' TH VERTEX NO LONGER WORST 
C BY DISPLACING TOWARD OVER-REFLECTING 

C THRU THE BEST VERTEX. 

IER = 2 

IF (NPR .LE. 0) GO TO 42 
WRITE (6,52) NT, J 

CALL BOUT (NT , NPT , NFE , NCE , NV , NVT , V ,K , FUN , CEN , J ) 
GO TO 42 

33 NT = NT+ 1 
GO TO 20 

C 

C SUCCESS: WE HAVE A REPLACEMENT FOR VERTEX J. 

34 FUN(J) = FUNTRY 
FUNOLD = FUNTRY 
NPT = NPT+1 



C 

C EVERY 100 ’TH PERMISSIBLE TRIAL, RECOMPUTE 
C CENTRIOD SUMMATION TO AVOID CREEPING ERROR. 

IF (MOD (NPT , 100 ) . NE . 0 ) GO TO 37 
C 

DO 36 1=1, NV 
SUM ( I ) = 0. 

C 

DO 35 N= 1 ,K 

35 SUM(I) = SUM(I ) +V (I , N ) 

C 

CEN ( I ) = SUM (I ) /FK 

36 CONTINUE 



LC = 0 
GO TO 39 
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c 

37 DO 38 1=1, NV 

38 SUM(I) = SUM( I ) +V( I , J ) 

C 

LC = J 
C 

39 IF (NPR.LE.O) GO TO 40 

IF (MOD(NPT,NPR) .NE.O) GO TO 40 
C 

CALL BOUT (NT , NPT , NFE , NCE , NV , NVT , V ,K , FUN , CEN , LC ) 
C 

C HAS THE MAX. NUMBER OF TRIALS BEEN REACHED 
C WITHOUT CONVERGENCE? 

C IF NOT, GO TO NEW TRIAL. 

40 IF (NT.GE.NTA) GO TO 41 
C 

C NEXT- TO -WORST VERTEX NOW BECOMES WORST. 

J = JN 
GO TO 17 

41 IER = 3 

IF (NPR.GT.O) WRITE (6,54) 

C 

C COLLECTOR POINT FOR ALL ENDINGS. 

C 1) CANNOT DEVELOP FEASIBLE VERTEX. IER = 1 

C 2) CANNOT DEVELOP A NO-LONGER- WORST VERTEX. IER = 2 

C 3) FUNCTION VALUE UNCHANGED FOR K TRIALS. IER = 0 

C 4) LIMIT ON TRIALS REACHED. IER = 3 

C 5) CANNOT FIND FEASIBLE VERTEX AT START. IER = -1 

42 CONTINUE 
C 

C FIND BEST VERTEX. 

CALL FBV (K , FUN ,M) 

IF (IER.GE.3) GO TO 44 
C 

C RESTART IF THIS SOLUTION IS SIGNIFICANTLY BETTER 
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C THAN THE PREVIOUS, OR IF THIS IS THE FIRST TRY. 

IF (NPR.LE.O) GO TO 43 
WRITE (6,55) (M , YMN , FUN (M) ) 

43 IF (FUN(M) .GE. YMN) GO TO 47 

IF (ABS(FUN(M)-YMN) . LE . AMAXl ( EP , EP*YMN ) ) GO TO 47 
C 

C GIVE IT ANOTHER TRY UNLESS LIMIT ON TRIALS REACHED. 

44 YMN = FUN (M) 

FUN ( I ) = FUN (M) 

C 

DO 45 1=1, NV 
CEN(I) = V(I,M) 

SUM (I ) = V ( I , M ) 

45 V(I , 1) = V(I,M) 

C 

DO 46 1=1, NVT 

46 XS(I) = V ( I ,M) 

C 

IF (IER.LT.3) GO TO 6 

47 IF (NPR.LE.O) GO TO 48 

CALL BOUT (NT , NPT , NFE , NCE , NV , NVT , V , K , FUN ,V(1,M) ,-l) 
WRITE (6,56) FUN(M) 

48 RETURN 
C 

49 FORMAT ( 50H0INDEX AND DIRECTION OF 
10UTLYING VARIABLE AT STARTI5 ) 

50 FORMAT (50H0IMPLICIT CONSTRAINT 
1VIOLATED AT START . DEAD END . ) 

51 FORMAT ('OCANNOT FIND FEASIBLE ’, 14 ,’ TH 
1VERTEX OR CENTROID AT START.’) 

52 FORMAT (10H0AT TRIAL I4,54H CANNOT FIND 
1FEASIBLE VERTEX WHICH IS NO LONGER 
1WORST, 14, 15X, ’RESTART FROM BEST VERTEX.’) 

53 FORMAT (40H0FUNCTION HAS BEEN ALMOST 
1UNCHANGED FOR 15, 7H TRIALS) 
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54 FORMAT (27H0LIMIT ON TRIALS EXCEEDED. ) 

55 FORMAT ( ' OBEST VERTEX IS NO. ',13,' 

1 OLD MIN WAS ' E15 . 7 ' , NEW MIN IS ’,E15.7) 

56 FORMAT ( ' OMIN OBJECTIVE FUNCTION IS ’,E15.7) 

END 

SUBROUTINE FBV (K , FUN ,M) 

DIMENSION FUN (50) 

M = 1 
C 

DO 1 1=2, K 

IF (FUN(M) .LE.FUN(I) ) GO TO 1 
M = I 

1 CONTINUE 
C 

RETURN 

END 

SUBROUTINE BOUT (NT , NPT , NFE , NCE , NV , NVT , V ,K , FN , C , IK ) 
DIMENSION V(50, 50) , FN(50), C(25) 

WRITE (6,4) NT, NPT, NFE, NCE 
C 

DO 1 1=1, K 

WRITE (6,5) FN(I) , ( V ( J , I ) , J = 1 , NV ) 

IF (NVT.LE.NV) GO TO 1 
NVP = NV+1 

WRITE (6,6) (V(J,I) , J=NVP,NVT) 

1 CONTINUE 
C 

IF (IK.NE.O) GO TO 2 
C 

WRITE (6,7) (C(I) ,I=1,NV) 

RETURN 

2 IF (IK.GE.O) GO TO 3 
WRITE (6,8) (C(I) ,I=1,NV) 

RETURN 

3 WRITE (6,9) IK, (C(I) ,I=1,NV) 
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RETURN 



4 FORMAT ('ONO. TOTAL TRIALS = ',I5,4X, 

1 ' NO . FEASIBLE TRIALS = ',I5,4X,'NO. FUNCTION 
1 EVALUATIONS = ',I5,4X,'N0. CONSTRAINT EVALUATI 
IONS = ' , 15/ ' 0 FUNCTION VALUE' ,6X, 'INDEPENDENT 

1VARIABLES/ DEPENDENT OR IMPLICIT CONSTRAINTS') 

5 FORMAT (1H , E18 . 7 , 2X , 7E14 . 7 / ( 21X , 7E14 . 7 ) ) 

6 FORMAT (21X.7E14.7) 

7 FORMAT (10H0CENTROID 11X , 7E14 . 7 / (21X , 7E14 . 7 ) ) 

8 FORMAT (’0 BEST VERTEX' , 7X, 7E14 . 7/ (21X, 7E14 . 7 ) ) 

9 FORMAT ( ' OCENTROID LESS VX' , 12 , 2X , 7E14 . 7 / 

1 (21X, 7E14 . 7 ) ) 

END 

FUNCTION FE (X ) 

DIMENSION X( 5 ) 

COMMON TDIFF 
CALL PLANT (X) 

FE=TDIFF 

RETURN 

END 

FUNCTION KE(X) 

DIMENSION X(5 ) 

KE = 0 

RETURN 

END 

//GO.SYSIN DD * 

/* 
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APPENDIX B 



CONTROLLER DESIGN 



A. CONTROLLER "C" 




Figure B.l Block Diagram of Controller C 
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Figure B.l corresponds to controller 'C' which has the 
form : 



Kjd + TiS) (1+T 4 S) 
(1+T 2 S) (I+T3S) 

Verifying that equation B.l 
' C ' we have : 



(B.l) 

corresponds to controller 



A= 



* 1*1 
i + t 2 s 



+ ^ i T i 2 



(B.2) 



^2_- l DX 2 =SX 2 

dx 2 s 



(B.3) 



X 2 _1 x = ^1 

?7 = i+t 2 s 2 'i+t 2 s 

1 (B . 4 ) 

Substituting equations B.3 , B.4 into equation B.2 we 

have : 



X 1 K 1 + K 1 T 1 X 1 S 
1 + T 2 S 




We also have: 



(B.5) 



X u = 



1 +ToS 



+ T4DX3 



(B . 6 ) 



^3 __1_ DX 3 = X 3 S 

dx 3 s 



(B • 7 ) 
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(B . 8 ) 



x 3 1 
A = 1+T 3 S 



•3 = 



A 



1 + TjS 



Substituting equations B.7 , B.8 into equation B.6 we 

have : 






A + T 4 AS 
1 + T 3 S 3 +T 3 S 



(B.9) 



Now substituting equation B.5 into equation B.9 we 



have : 



X = 



XlKi 



X1K1T1S 



(1+T 2 S) ( I+T 3 S) (1+T 2 S)(1+To 
Xi^TijS X 1 K 1 T 1 T,,s 2 



S) 



(1+T 2 S)(1+T 3 S) ( 1+T 2 S) ( i+t 3 s) 



(B.10) 



Finally rearranging terms equation B.10 can be 
written: 



X 4 _ K 1 (l^T 1 S-H t ,S-^T 1 T 4 S 2 ) _ K 1 (1 + T 1 S)(1 + T 4 S) 
X 1 (1+T 2 S) (I+T 3 S) a+T 2 S)(l+T 3 S) 



(B.ll) 



B. CONTROLLER "B” 

Figure B.2 corresponds to controller 'B' which has the 
form: 
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Figure B.2 Block Diagram of Controller B 



KiU + TiS) 
(1+T 2 S)(1+T 3 S) 



(B . 12 ) 



Verifying that equation B . 12 corresponds to Controller 
' B ' we have : 



_ K 1 X 1 

T+T 9 S + X 1^1^^2 

(B . 13 ) 
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1 



dx 2 = X 2 S 



(B.14) 



^2 

dx 2 ‘ S 



* 2 _ ! X 1 
X 1 't 2 s+1 X2 ' 1+T 2 S 



(B.15) 



Substituting equations B.14 , B.15 into equation B . 13 

we have : 



K 1 X 1 K iT 1 X 1 S 
A = + 

i+t 2 s i+t 2 s 



(B. 16 ) 



In a similar way we also can derive equation B.17 



X3 = 



A 

i + t 3 s 



(B.17) 



Substituting equation B.16 into equation B.17 we have: 



X3 = 



* 1 X 1 + 
(1+T 2 S) (1+T 3 S) 



KiTiXiS 

( 1 +T 2 S) (I+T 3 S) 



(B. 18) 



Finally rearranging terms equation B.18 becomes: 



X 3 Kid + TjS) 
x7 = (l+T 2 S)(l+T 3 S) 



(B.19) 
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C. CONTROLLER "A” 




Figure B.3 Block Diagram of Controller A 



Figure B.3 corresponds to controller 'A' which has the 
f orm : 

Kxd+TjS) 

l + T-S 

2 (B . 20 ) 
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Verifying that equation B.20 corresponds to Controller 
'A' we have: 



, X 1 K 1 
i+t 2 s 



K 1 t 1 dx 2 



(B . 21 ) 



But since we know that: 




dx 2 =x 2 s 



(B . 22 ) 



*2. 1 y , X 1 

Xj^ HT 2 S 1*T 2 S (b.23) 



Substituting equation B.22 , B.23 into equation B.21 

we have : 



*3 = 



X 1 K 1 K 1 T 1 X 1 S 
+ 

i+t 2 s i+t 2 s 



(B. 24) 



Finally rearranging terms equation B . 24 becomes: 



X 3 K 1 (1+T 1 S) 

1+T 2 S (B . 25 ) 
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D. CODING OF THE EQUATIONS 



For Controller "C” 



Integration 



YAWE = YAWC - YAW 

YAWE = YAW - YAWC 

CYAWE - X 2 ) 

DX 9 = — 

T 2 

A = K x (X 2 + T 1 DX 2 ) 
A-X o 

dx 3 = 

t 3 

D = X 3 + T 4 DX 3 



X 2 = x 2 + DX 2 - delt 
X 3 = x 3 + DX 3 - delt 



For Controller "B" Integration 



YAWE = YAW 
YAWE 


- YAWC 

- X~ 






DX2 = 


2 


X 2 = x 2 + dx 2 


• DELT 


T 2 

A = K 3 (X 2 + 


TiDX 2 ) 


x 3 = x 3 + DX 3 


• DELT 



For Controller "A" Integration 

YAWE = YAW - YAWC 
YAWE - X. 

DX 0 = 2 - X 9 = X 9 + DX 9 • DELT 

2 T 2 1 1 l 

D - K^(X 2 + T 3 DX 2 ) 
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For all of the above cases the equations include the 
error detector and the controller which are indicated in 
Figure B.4 




Figure B.4 General Scheme of Control 
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APPENDIX C 



RESPONSE OF THE SYSTEM FOR REGULAR SEAS 

//PROGRA JOB (???? ,0356) , 'RESEARCH' ,CLASS=A 
//"MAIN 0RG=NPGVM1. ????P 

// EXEC FORTXCG , PARM . FORT= ' OPT (2 ) ' , IMSL=DP , REGION= 1024K 
//FORT. SYS IN DD * 

C 

C 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS 
C HAVE BEEN OBTAINED CHANGE XS(*) TO X(*) AND DELETE 
c XU ( * ) , AND XL("). 

COMMON J 
DIMENSION X(5 ) 

X(l)=l. 8287125 
X(2)=l. 1652012 
X(3)=10. 5659571 
X(4)=ll. 7124157 
X(5)=20. 5683947 
C CALL PLANT (X) 

C IF ONLY SIMULATION IS WANTED 
CALL PLANT (X) 

WRITE (6,25) 

25 FORMAT (IX,' OPTIMAL GAINS ',/ ) 

DO 30 1=1,5 

30 WRITE (6,40)I,X(I) 

40 FORMAT (IX, 'X(' ,12, ' )=' , F14 . 7 ) 

WRITE (6, 50) J 

50 FORMAT ( IX, ’J = ' ,E15.10) 

STOP 

END 

SUBROUTINE PLANT (XX) 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
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COMMON TDIFF 

REAL, *8 L,L2 ,L3 , L4 ,L5 ,L6 

REAL* 8 X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL* 8 TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL* 8 YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
REAL* 8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL* 8 RHO , IZ , FX , FY , MZ , XP , MASS , DELT , MZI , RXI , WA , WE 
REAL* 8 DYAWE , YAWE , YAWC , I SE , I SR , LAMDA , D , RYR , RYI , MZR 
REAL* 8 K1 ,Tl,T2 ,D,X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR 
REAL* 8 T3 ,T4 ,X3 , DX3 , X4 
DIMENSION XX (5) 

C 

C CLOSE LOOP ANALYSIS WITH FILTER 

C 

C INITIAL CONDITIONS FOR INTEGRATION 

C SIMULATION END TIME IN SECONDS 
ETIME=600 . 

TIME=0 . 0 
I COUNT = 1 . 0 

C INITIALIZE THE COST FUNCTION 
ISE=0 . 0 
ISR=0 . 0 
TDIFF=0 . 0 
LAMDA =8 . 128 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
K1=XX ( 1 ) 

T1=XX(2) 

T2=XX (3 ) 

T3=XX (4 ) 

T4=XX(5 ) 

C WRITE(6 , 1010) K1,T1,T2 

C1010 FORMAT ( IX , ' K1 =',F15.7,'T1 =',F15.7,’T2 =',F15.7) 

C X, XDOT, Y, YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0 . 0 
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XDOT=0 . 0 
YDOT=0 . 0 



C U , UDOT , V , VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UD0T=0 . 0 
VDOT=0 . 0 
YAW=0 . 0 
R=0 . 0 
RDOT=0 . 0 

C ORDERED SPEED IN FEET/ SEC 
C 38.82 FT/ SEC=23 KNOTS 
UC=38 .82 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 

C D = RUDDER ANGLE 
D=0 . 0 
L=880 . 5 
L2=L**2 
L3=L*L*L 
L4=L*L3 
L5=L*L4 
L6=L*L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
C MOMENTS IN Z 
FX=0 . 

FY=0 . 

MZ = 0 . 

C RXR= -0 . I5744D+05 

C RXI=-0 .19950D+06 

C RYR=0 .52365D+04 

C RYI=0. 18699D+06 

C MZR= -0 .29870D+08 

C MZI=-0 .35751D+07 

C RXR= -0 . 50230D+04 
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c 



c 



c 



c 



c 



RXI=0. 127 12D+05 
RYR=0 . 35290D+04 
RYI= -0 . 31909D+05 
MZR=0 . 38826D+07 
MZI=-0 . 643 13D+07 



RXR=0 . 28540D+04 
RXI = - 0 . 99574D+04 
RYR=-0 . 85441D+04 
RYI=0 .39595D+05 
MZR=-0 . 13014D+08 
MZI=0. 11348D+08 
C RXR= -0 .75642D+04 

C RXI=0 . 83497D+04 

C RYR=0 .23379D+05 

C RYI = - 0 . 8 1502D+ 05 

C MZR=0 . 28622D+07 

C MZI=-0 . 19388 D+ 08 

C RXR=-0 .37916D+04 

C RXI=0 . 16381D+04 

C RYR=-0 . 76647D+05 

C RYI = -0 . 37 685D+05 

C MZR=-0 .83915D+07 

C MZI=-0 . 53176D+07 

RX=DSQRT ( RXR**2 +RXI**2 ) 

RY =DSQRT ( RYR**2 + RYI**2 ) 

RZ=DSQRT ( MZR**2 +MZ 1**2 ) 

TX=DATAN (RXI/RXR) 

T Y = DAT AN ( RY I / R YR ) 

TZ=DATAN (MZI /MZR) 

C SIGNIFICANT WAVE HEIGHT; SEA STATE 1-5,2-10,3-15, 

C 4-17.5,5-22.5 6-27,7-35,8-42,9-60 
WA= 17 . 5 

C ENCOUNTER FREQUENCY .1,. 2, .3, .4, .5, .6, .75, 1.0, 1.5, 2. 5 
WE=0 . 6 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE 
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c AS PARAMETERS 
RH0=1. 9876 

MASS= ( . 0044 )* ( . 5 "RHO“L3 ) 

IZ=(0. 00028 )*(.5*RHO*L5) 

YAWE = 0 . 0 
X2=0 . 0 
DX2=0 . 0 
X3 = 0 . 0 
DX3=0 .0 
X4=0 . 0 

200 CONTINUE 

S=DSQRT (u**2+v**2 ) 

C INPUT YAW COMMAND 
YAWC=0 .0 

IF (TIME. GE. 0.0) YAWC=0.0 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL- YAW ORDERED) 
C ( COMPENSATOR FILTER ) 

YAWE= YAW - YAWC 
DX2= ( YAWE-X2 ) /T2 
X4=K1*(T1*DX2+X2) 

DX3= (X4-X3 )/T4 
D=(T3*DX3+X4) 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

C XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER 
C FREQUENCY 

XUDOT= ( - . 000 1 )* ( . 5*RHO*L3 ) 

XU= (-0 .0253)*( .5*RHO*L2*S) 

XUU= ( -0 . 0003 )*( . 5*RHO*L2) 

XVR= ( 0 . 00 3 9 )* ( . 5*RHO*L3 ) 

XVV= ( - . 00 12 ) * ( . 5*RHO*L2 ) 

XDD = ( - 0 . 0 0 0 5 ) * ( . 5 "RHO *L2 * S ** 2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

C YV= ( -0 . 00758 )*( . 5*RHO*L2*S ) 

YR= (0 . 0023 ) * ( . 5*RHO*L3*S ) 



99 



YD= (0 . 00145 )*( . 5*RHO*L2*S**2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L 3 / S ) 

YVRR= ( - 0 . 008 )* ( . 5 "RH0“L4 / S ) 

YVVV= ( - 0 . 0 3 ) * ( . 5 "RHO"L2 / S ) 

YRRR= ( 0 . 0 0 3 ) ' * ( . 5 *RHO*L5 / S ) 

YDDD= ( - 0 . 0005 )* ( . 5*RH0*L2*S**2 ) 

C YUDOT IS THE ADDED MASS TERM WHICH MUST BE 
C CHANGED FOR DIFFERENT ENCOUNTER ANGLE , SPEED , 

C ENCOUNTER FREQUENCY 
C YVDOT= ( - 0 . 003 9 ) * ( . 5*RHO*L3 ) 

C SPEED=23 KNOTS, ENCOUNTER ANGLE = 30 , ENCOUNTER 
c FREQUENCY =0.4 
C YVD0T=-0. 30908D+07 

C YV= - 0 . 8 127 1D+04 

C YVD0T=-0 . 36185D+07 

C YV=-0 . 24757D+06 

YVDOT = -0.32890D + 07 
YV=-0 . 11775D+07 
C YVDOT= -0.23038D+07 

C YV=-0 . 18267D+07 

C YVDOT= -0.59800D+06 

C YV=-0 . 13260D+07 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV= ( - 0 . 00213 ) * ( . 5*RHO*L3*S ) 

C NR= ( - 0 . 00105 )* ( . 5*RHO*L4*S ) 

ND= ( - 0 . 000 7 ) * ( . 5 *RHO *L 3 * S * * 2 ) 

NVVR= ( - 0 . 0 15 ) * ( . 5 *RH0*L4 / S ) 

NVRR= ( - 0 . 0 0 8 ) * ( . 5 *RHO*L5 / S ) 

NVVV= ( 0 . 0 1 ) * ( . 5 "RHO"L3 / S ) 

NRRR= ( - 0 . 006 )* ( . 5-'RHO-'L6 / S ) 

NDDD= (0 . 0001)*( . 5 *RHO *L 3 * S * * 2 ) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER 
C FREQUENCY 

C NRDOT= (-0 . 00027 )*( .5*RHO*L5) 
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C SPEED=23 KNOTS, ENCOUNTER ANGLE = 30 , ENCOUNTER 
C FREQUENCY =0.4 
C NRDOT=-0 .2625 ID +12 

C NR= -0 .53637D+09 

C NRDOT= -0.20125D+12 

C NR= -0.94970D+10 

NRDOT=-0. 18671D+ 12 
NR= -0 .46860D+11 
C NRDOT=-0 . 145 18 D+ 12 

C NR= -0.87538D+11 

C NRDOT=-0 .3726 ID +11 

C NR= -0 . 69856D+ 11 

C REGULAR WAVE SEA STATE 

FX=WA*RX*DCOS (WE*TIME+TX ) 

FY = WA"RY“DCO S ( WE "TIME + TY ) 

MZ=WA*RZ*DCOS (WE*TIME+TZ ) 

C U ACTUAL SPEED 
C UC COMMANDED SPEED 
C XP = PROPELLER THRUST 
X P = - XUU *U C * * 2 
C EQUATIONS OF MOTION 

C UDOT= ( (XVR + MASS ) *V*R + XUU*U**2 + XVV*V**2 

C 1 + XDD"D"D + FX + XP ) / (MASS-XUDOT ) 
VDOT=(YV*V + (YR-MASS"U)"R + YD*D 
1 + YVVR*V**2*R + YVRR*V*R**2 
1 + YVW*V**3 + YRRR"R" "3 + YDDD*D**3 

1 + FY )/ (MASS-YVDOT) 

RDOT= ( NV"V + NR*R + ND"D + NWR*V**2*R 
1 + NVRR*V*R**2 

1 + NVW*V**3 + NRRR*R**3 + NDDD*D**3 
1 + MZ ) / (iz-NRDOT ) 

C WHEN TO PRINTOUT 

IF (ICOUNT.EQ. 2) GO TO 50 
GO TO 300 

C CONVERT RADIANS TO DEGREES 
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50 YAWDEG= YAW*57.296 
RDEG=R“57 .296 
RDDEG=RDOT*57 .296 
DDEG=D*57 .296 
YAWC=YAWC"57 .296 
WRITE (6,100) TIME , YAWDEG 

C 1 , UC , U , UDOT , V , VDOT , YAWC , YAWDEG , RDEG , RDDEG , DDEG 

100 FORMAT ( IX , F 12 . 8 , IX , F 12 . 8 ) 

C 1' FT XDOT= ' ,F8 .4, ' FT/ SEC Y=',F8.2,' FT 
C 1 YDOT = ’ , F8 . 4 , ' FT/ SEC 

C 1 ' , / , 2X , ' UC= ' , F8 . 4 , ' FT/SEC U=',F8.4,' 

C 1 FT/SEC UDOT= ' , F10 . 6 , 

C 1 ’ FT/ SEC ''*2 V= ’ , F8 . 4 , ’ FT/ SEC 

C 1 VDOT= ’ ,F10. 6 , ’ FT/SEC**2’ 

C 1 , / , 2X , ' YAWC= ' , F8 . 4 , ' DEG YAW= ' 

C 1 , F15 . 7 , ' DEG YAW RATE= ' , F15 . 7 , ' 

C 1 DEG/SEC YAW ACCEL= ’ 

C 1 , F15 . 7 , ' DEG/SEC**2’ ,/ ,2X, ' 

C 1 RUDDER = ' , F15 . 7 , ' DEG ') 

C WRITE (6,101) TIME, DDEG 

C101 FORMAT ( IX , F12 . 8 , IX , F12 . 8 ) 

I COUNT =1 

C TEST IF WANT TO STOP 
300 IF ( TIME . GE . ETIME ) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELT= 1 . 0 
C INTEGRATION 

U=U+UDOT*DELT 
V=V+VDOT*DELT 
R=R+RDOT*DELT 
YAW= YAW + R*DELT 
X2=X2+DX2*DELT 
X 3 = X 3 + DX 3 *DELT 

C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT=U*DCOS (YAW) -V"DSIN (YAW) 
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YDOT=U*DSIN (YAW) +V*DCOS ( YAW) 

X=X+XDOT*DELT 
Y = Y + YDO T *DELT 
TIME=TIME+DELT 
ICOUNT=ICOUNT+ 1 
ISE=ISE + LAMDA*YAWE**2 
ISR=ISR + D**2 
GO TO 200 

C J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 

WRITE (6, 500) ISE , ISR, TDIFF ,K1,T1,T2,T3 ,T4 
500 FORMAT (’ ' , IX , ' ISE= ' , F15 . 7 , ' ISR= ' , F15 . 7 , ' 

1 TOTAL= ' , F15 . 7 , 2x , 

1 'Kl=' ,F15.7,2X, 'Tl=' ,F15.7,2X, 'T2=' , F15 . 7 , 2X , 
1 'T3=' ,F15.7,2x, ' T4= ' ,F15.7) 

RETURN 

END 

//GO. SYS IN DD * 

/* 
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APPENDIX D 



DETERMINATION OF OPTIMAL CONTROLLER PARAMETERS FOR IRREGULAR 

SEAS 

/ /TRIALl JOB (1707 ,0356) , 'RESEARCH' ,CLASS=C 

//"MAIN 0RG=NPGVM1. 1707P 

// EXEC FORTXCG , PARM . FORT= ' OPT (2 ) ' , IMSL=DP , REGION= 1024K 

/ / FORT . SYSIN DD * 

C THIS PROGRAM WILL OBTAIN THE CONTROLLER OPTIMAL 

C GAINS. IT IS REFERENCED IN CHAPTER 5. 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS 

C HAVE BEEN OBTAINED CHANGE XS(*) TO X ( " ) AND 

C DELETE XU ( * ) , AND XL ( * ) . 

DIMENSION XS(5) ,XU(5) ,XL(5) 

XS(1)=0. 655751 
XS (2)=80 . 5483 
XS(3)=10. 74847 
XS (4 ) = 12 . 9 
XS (5 )=45 . 09 

C XS(I) IS THE STARTING GUESS 

C XL (I) IS THE LOWER LIMIT FOR THE I ' TH VARIABLE 

C XU (I) IS THE UPPER LIMIT FOR THE I ’ TH VARIABLE 
XL ( 1 ) = 0 . 1 
XU ( 1 ) = 2 . 5 
XL(2 ) =40 . 0 
XU(2 ) = 100 . 0 
XL ( 3 ) = 0 . 1 
XU ( 3 ) =20 . 0 
XL ( 4 ) = 5 . 0 
XU ( 4 ) = 8 0 . 0 
XL ( 5 ) = 6 0 . 0 
XU(5 ) = 150 . 0 

C A DESCRIPTION OF THE FOLLOWING PARAMETERS 
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C IS DISCUSSED IN BOXPLX 
R=9. / 13 . 

NTA= 1000 
NPR= 100 
NAV=0 
NV=5 
IP = 0 

C THE FOLLOWING STATEMENT MUST BE CHANGED TO 
C CALL PLANT (X) 

C IF ONLY SIMULATION IS WANTED 

CALL BOXPLX ( NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER ) 
WRITE (6,25) 

25 FORMAT (IX,' OPTIMAL GAINS ',/ ) 

DO 30 1=1,5 

30 WRITE (6,40)I,XS(I) 

40 FORMAT (IX, 'X(' ,12, ' )=' ,F14.7) 

STOP 

END 

SUBROUTINE PLANT (XX) 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
COMMON TDIFF 
REAL* 8 L,L2,L3,L4,L5,L6 

REAL*8 X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 

REAL*8 TIME , ETIME ,XUDOT , XUU , XVR , XVV , XDD 

REAL *8 YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 

REAL*8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL*8 RHO , IZ , FX , FY , MZ , XP , MASS , DELT 

REAL* 8 DYAWE , YAWE , YAWC , I S E , I SR , LAMDA , D 

REAL* 8 Kl,Tl,T2,T3,T4,D,X2, DX2 , X3 , DX3 , X4 , CH ( 1 1 ) , S 

DIMENSION XX (5 ) 

C 

C CLOSE LOOP ANALYSIS WITH FILTER 
C 

C INITIAL CONDITIONS FOR INTEGRATION 
C SIMULATION END TIME IN SECONDS 
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ETIME=600 . 

TIME=0 . 0 
ICOUNT= 1 

C INITIALIZE THE COST FUNCTION 
ISE=0 . 0 
ISR=0 . 0 
TDIFF=0 . 0 
LAMDA=8. 128 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
K1=XX(1) 

T1=XX(2) 

T2=XX ( 3 ) 

T3=XX (4 ) 

T4=XX(5 ) 

C WRITE(6 , 1010) K1,T1,T2 

C1010 FORMAT ( IX , ' K1 =',F15.7,' Tl =’,F15.7,' 

C 1 T2= ’ , F15 . 7 ) 

C X , XDOT , Y , YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0 . 0 • 

XD0T=0 . 0 
YD0T=0 . 0 

C U ,UDOT ,V , VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UD0T=0 . 0 
VDOT=0 . 0 
YAW=0 . 0 
R=0 . 0 
RD0T=0 . 0 

C ORDERED SPEED IN FEET/ SEC 

C 38.82 FT / SEC=23 KNOTS 
UC=38 . 82 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 

C D = RUDDER ANGLE 
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D=0 . 0 

L=880 . 5 

L2=L**2 

L3=L*L*L 

L4=L*L3 

L5=L*L4 

L6=L*L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
C MOMENTS IN Z 
FX=0 . 

FY=0 . 

MZ = 0 . 

C I SEA IS A SWITCH ; ISEA=0 (CALM WATER) ISEA=1 (SEA STATE) 
ISEA= 1 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE AS 
C PARAMETERS 

RHO = 1.9876 

MASS= ( . 0044)* ( . 5*RHO*L3 ) 

IZ=(0. 00028 )*(.5*RHO*L5) 

YAWE=0 . 0 
X2 = 0 . 0 
DX2=0 . 0 
X3=0 . 0 
DX3=0 . 0 
X4 = 0 . 0 

200 CONTINUE 

S=DSQRT (U**2 +V**2 ) 

C INPUT YAW COMMAND 
YAWC=0 . 0 

IF (TIME. GE. 0.0) YAWC=0.0 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL- YAW ORDERED) 

C ( CONTROLLER FILTER ) 

YAWE=YAW - YAWC 
DX2= ( YAWE-X2 ) /T2 
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X4 = K 1* ( T 1*DX2 + X2 ) 

DX3= (X4-X3 ) /T4 
D= (T3"DX3 + X4) 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

C XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER 
C FREQUENCY 

XUDOT = ( - .0001)*( .5*RHO*L3) 

XU= ( - 0 . 0253 )'*'(. 5*RHO*L2*S ) 

XUU= ( - 0 . 0003 ) * ( . 5*RHO*L2 ) 

XVR= ( 0 . 0 0 3 9 ) * ( . 5 *RHO*L3 ) 

XVV= ( - . 0012 )* ( . 5*RHO*L2 ) 

XDD= ( - 0 . 0005 ) * ( . 5 * RH O * L 2 * S * * 2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

YV= ( - 0 . 00758 ) * ( . 5*RHO*L2*S ) 

YR= (0 . 0023 )*( . 5*RHO*L3*S ) 

YD= (0 . 00145 )*( . 5 *RHO*L 2 * S * * 2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= ( - 0 . 0 0 8 ) * ( . 5 *RHO*L4 / S ) 

YVVV= ( - 0 . 03 ) * ( . 5*RHO*L2 / S ) 

YRRR= ( 0 . 0 0 3 ) * ( . 5 *RHO*L5 / S ) 

YDDD= ( - 0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 

C YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER 
C FREQUENCY 

C YVDOT= ( -0 . 003 9 )* ( . 5*RHO*L3 ) 

C SPEED=23 KNOTS, ENCOUNTER ANGLE = 

C ENCOUNTER FREQUENCY = . 75 
YVDOT= -2304300.0 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV= ( - 0 . 002 13 ) * ( . 5*RHO*L3*S ) 

NR= ( - 0 . 00 105 ) * ( . 5*RHO*L4*S ) 

ND= (-0 . 0007 )* ( . 5*RHO*L3*S**2) 

NVVR= ( - 0 . 0 15 ) * ( . 5*RHO*L4 / S ) 

NVRR= ( - 0 . 008 ) * ( . 5*RHO*L5 / S ) 
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c 

c 

c 

c 

c 

c 

c 



c 

c 



c 

c 

c 

c 

c 

c 



NVVV= ( 0 . 0 1 ) * ( . 5*RHO*L3 / S ) 

NRRR= ( - 0 . 00 6 ) * ( . 5 *RH0 *L6 / S ) 

NDDD= ( 0 . 000 1 )* ( . 5 * RH 0 * L 3 * S * * 2 ) 

NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE 
CHANGED FOR DIFFERENT ENCOUNTER ANGLE , SPEED 
ENCOUNTER FREQUENCY 

NRDOT= ( -0 . 00027 )* ( . 5*RHO*L5 ) 

SPEED=23 KNOTS, ENCOUNTER ANGLE = 

ENCOUNTER FREQUENCY =.75 
NRDOT= - 1 . 45 18E+ 11 
SETS SEA STATE TO ZERO 

IF (ISEA.EQ. 1) GO TO 30 
FX=0 . 

FY=0 . 



MZ = 0 . 

GO TO 35 

UNIT 12 HAS THE SEA STATE DATA NAMED CH 
IT MUST BE SYNCHRONIZED BY TIME 
30 READ (12) CH 
FX=CH( 3 ) 

FY=CH(4) 

MZ=CH( 8 ) 

35 CONTINUE 
U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP= - XUU*UC**2 



EQUATIONS OF MOTION 

UDOT= ( (XVR + MASS )*V*R + XUU*U**2 + XVV*V 
1 + XDD*D*D + FX + XP )/ (MASS-XUDOT) 

VDOT= (YV*V + ( YR- MAS S "U ) *R + YD*D + YVVR*V 



1 + YVRR*V*R**2 

1 + YVW*V**3 + YRRR*R**3 

1 + FY )/ (MASS-YVDOT) 

RDOT=( NV*V + NR"R + ND-'D 



YDDD*D**3 

NVVR*V**2*R 



2 

2*R 
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1 + NVRR*V*R**2 

1 + NVVV*V**3 + NRRR*R**3 + NDDD*D**3 
1 + MZ ) / (IZ-NRDOT ) 

C WHEN TO PRINTOUT 

IF (ICOUNT.EQ.il) GO TO 50 
GO TO 300 

C CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW" 57. 29 6 
RDEG=R*57 .296 
RDDEG=RDOT*57 .296 
DDEG=D*57 .296 
YAWC = YAWC*5 7.296 

C WRITE (6,100) TIME , XP , X , XDOT , Y , YDOT 

C 1 , UC , U , UDOT , V , VDOT , YAWC , YAWDEG , RDEG , RDDEG , DDEG 

100 FORMAT ( IX , ' TIME= ' , F8 . 3 , ' SEC XP=’,F10.2,’ LBF 
1 X=',F8.2, 

1' FT XDOT= ’ ,F8 .4, ' FT/ SEC Y=’,F8.2,’ FT YDOT= ' 
1 , F8 . 4 , ' FT/SEC 

r,/,2X,’ UC= ' , F8 . 4 , ' FT/SEC U=’,F8.4,’ FT/SEC 
1 UDOT = '* , F10 . 6 , 

1 ' FT / SEC**2 V= ' , F8 . 4 , ' FT/ SEC VDOT= ' , F10 . 6 , 

1 ’ FT / SEC" *2 ' , / , 2X , ' YAWC= ' , F8 . 4 , ' DEG YAW= ' 

1 , F15 . 7 , ' DEG YAW RATE= ’ ,F15 .7 , ' DEG/ SEC 

1 YAW ACCEL= ’ 

1 , F15 . 7 , ' DEG / SEC**2 ’ , / , 2X , ' RUDDER = ' , F15 . 7 , ' 

1 DEG ' ,/) 

I COUNT =1 

C TEST IF WANT TO STOP 
300 IF (TIME . GE . ETIME ) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELT=1 . 0 
C INTEGRATION 



U=U+UDOT"DELT 

V=V+VDOT"DELT 

R=R+RDOT"'DELT 
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YAW=YAW+R*DELT 
X 2 = X2 + DX 2 *DELT 
X 3 = X 3 + DX 3 *DELT 

C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
C XDOT=U' v DCOS ( YAW) -V*DSIN ( YAW) 

C YDOT=U*DSIN(YAW)+V*DCOS(YAW) 

C X=X+XDOT*DELT 

C Y = Y+ YDOT *DELT 

TIME=TIME+DELT 
ICOUNT=ICOUNT+ 1 
ISE=ISE + LAMDA-'YAWE 2 
ISR=ISR + D**2 
GO TO 200 

C J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 

WRITE (6 ,500) TDIFF ,K1,T1,T2,T3,T4 
500 FORMAT(' ’ ,IX, 'TDIFF =’ ,F15.7, ’ K1=',F15.7,’ 

1 Tl = ' , F15 . 7 ,2X, 

1 ' T2= ' , F15 . 7 , 2X , ' T3= ' , F15 . 7 . 2X , ’T4=’ ,F15.7) 
REWIND 12 
RETURN 
END 

C BETWEEN LINE 249 (END) AND THE FOLLOWING LINE 

C (//GO.SYSIN DD *)WE HAVE TO INCLUDE BOXPLX . 

//GO.SYSIN DD * 

/* 

/ /GO . FT12F001 DD DISP=SHR ,DSN=MSS . S2160 . A213 
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APPENDIX E 



RESPONSE OF THE SYSTEM FOR IRREGULAR SEAS 

//PROGRA JOB (???? ,0356) , 'RESEARCH' ,CLASS=B 
//"MAIN 0RG=NPGVM1. ????P 
// EXEC FRTXCLGP , IMSL=DP , REGION^ 1024K 
//FORT. SYS IN DD * 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS 
C HAVE BEEN OBTAINED. 

DIMENSION XX ( 5 ) 

C OPTIMAL GAINS FOR CONTROLLER 
XX(1)=2. 45967680 
XX(2)=88. 2797241 
XX(3)=50. 5678864 
XX(4)=5. 27039050 
XX(5)=95. 3189392 

C THE SUBROUTINE PLANT SIMULATES THE SL-7 CONTAINERSHIP 
CALL PLANT (XX) 

WRITE (6, 25) 

25 FORMAT (IX, 'OPTIMAL GAINS' ,/) 

DO 30 1=1,5 

30 WRITE (6,40)1, XX ( I ) 

40 FORMAT (IX, 'XX( ' ,12, ' )= ' ,F14. 7) 

STOP 

END 

C 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
SUBROUTINE PLANT (XX) 

COMMON TDIFF 

REAL* 8 L,L2,L3 ,L4,L5 ,L6 

REAL*8 X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL* 8 TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL* 8 YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
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REAL* 8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL*8 RHO , IZ , FX , FY , MZ , XP , MASS , DELT 
REAL*8 DYAWE , YAWE , YAWC , I SE , I SR , LAMDA , D 
REAL* 8 K1,T1,T2,D,X2, DX2 , S , CH( 11 ) ,DX3,X3,X4 
DIMENSION XX(5 ) 

C 

C CLOSE LOOP ANALYSIS WITH FILTER 

C 

C INITIAL CONDITIONS FOR INTEGRATION 

C SIMULATION END TIME IN SECONDS 
ETIME=600 . 

TIME=0 . 0 
I COUNT =1 

C INITIALIZE THE COST FUNCTION 
ISE=0 . 0 
ISR=0 . 0 
TDIFF=0 . 0 
LAMDA= 4 . 2 

C GAIN COEFFICIENTS 
K1=XX( 1) 

Tl=XX(2 ) 

T2=XX( 3 ) 

T3=XX(4) 

T4=XX(5) 

C X , XDOT , Y , YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0 . 0 
XDOT=0 . 0 
YDOT=0 . 0 

C U , UDOT , V , VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT=0 . 0 
VDOT=0 . 0 
YAW= 0 . 0 
R=0 . 0 
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RD0T=0 . 0 
YAW=0 . 0 

C ORDERED SPEED IN FEET/ SEC 
C 54.01 FT/ SEC=32 KNOTS 
UC= 38 .81 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 

C D = RUDDER ANGLE 
D=0 . 0 
L=880 . 5 
L2=L**2 
L3=L*L*L 
L4=L*L3 
L5=L“L4 
L6=L*L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
C MOMENTS IN Z 
FX=0 . 

FY=0 . 

MZ = 0. 

C ISEA IS A SWITCH; ISEA=0(CAL WATER) ISEA= 1 ( SEA STATE) 
ISEA= 1 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE AS 
C PARAMETERS . 

RHO= 1.9876 

MASS= ( . 0044 )*( . 5*RHO*L3 ) 

IZ= (0 . 00028 )*( . 5*RHO*L5 ) 

YAWE = 0.0 
X2 = 0 . 0 
DX2=0 . 0 
X3=0 . 0 
DX3=0 . 0 
X4=0 . 0 

200 CONTINUE 
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S=DSQRT(U**2 + V**2) 

C INPUT YAW COMMAND 
YAWC=0 . 0 

IF (TIME . GE . 0 . 0 ) YAWC=0.0 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL- YAW ORDERED) 
C ( COMPENSATOR FILTER ) 

YAWE= YAW - YAWC 
DX2 = ( YAWE - X2 ) / T2 
X4=K1' V (T1*DX2 + X2 ) 

DX3= (X4-X3 ) /T4 
D= (T3' V DX3 + X4) 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

C 

C XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE AND SPEED. 

C XUDOT = ( - . 0001 )* ( . 5*RHO*L3 ) 

XUU= ( -0 . 0003 )*( . 5*RHO*L2 ) 

XVR= ( 0 . 0 0 3 9 ) * ( . 5 *RHO*L3 ) 

XVV= ( - . 00 12 )* ( . 5*RHO*L2 ) 

XDD= ( - 0 . 0005 )* ( . 5*RHO*L2*S**2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

YV= ( -0 . 00758 )*( . 5 "RHO“L2“ S ) 

YR= (0 . 0023 )*( . 5*RHO*L3*S ) 

YD= (0 . 00145 )*( . 5*RHO*L2*S**2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= ( - 0 . 0 08 ) * ( . 5 *RHO*L4 / S ) 

YVVV= ( - 0 . 0 3 ) * ( . 5 *RHO*L2 / S ) 

YRRR= ( 0 . 0 0 3 ) * ( . 5 *RHO*L5 / S ) 

YDDD- (-0.0005 )•’■'( . 5*RHO*L2*S**2 ) 

C YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE AND SPEED. 

C 

C YVDOT =(-0.0039) *(.5 *RHO*L3 ) 

YVDOT=-2304300 . 00 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
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NV= ( - 0 . 00213 ) * ( . 5*RHO*L3*S ) 

NR= ( - 0 . 0 0 1 0 5 ) * ( . 5 *RHO*L4*S ) 

ND= (-0 . 0007 )*( . 5*RHO*L3*S**2 ) 

NVVR= ( - 0 . 0 15 ) * ( . 5*RHO*L4 / S ) 

NVRR= ( - 0 . 0 0 8 ) * ( . 5 *RHO*L5 / S ) 

NVVV= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

NRRR= ( - 0 . 0 0 6 ) * ( . 5 *RHO*L6 / S ) 

NDDD= ( 0 . 000 1 ) * ( . 5 *RHO*L3 *S**2 ) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE AND SPEED. 

C 

C NRDOT= ( -0 . 00027 )*( . 5*RHO*L5 ) 

NRDOT=- 1 .5096E+11 
C SETS SEA STATE TO ZERO 

IF (I SEA. EQ. 1) GO TO 30 
FX=0 . 

FY = 0 . 

MZ = 0 . 

GO TO 35 

C UNIT 12 HAS THE SEA STATE DATA NAMED CH 
C IT MUST BE SYNCHRONIZED BY TIME 
30 READ (12) CH 

FX= CH(3 ) 

FY= CH(4) 

MZ= CH(8 ) 

35 CONTINUE 

C U ACTUAL SPEED 
C UC COMMANDED SPEED 
C XP = PROPELLER THRUST 
XP= -XUU*UC**2 
C EQUATIONS OF MOTION 

C UDOT= ( (XVR + MASS )*V*R + XUU*U**2 + XVV*V**2 
C 1 + XDD*D*D + FX + XP ) / (MASS-XUDOT ) 

VDOT= (YV*V + ( YR-MASS " S ) *R + YD*D + YWR*V**2*R 
1 +YVRR*V*R**2 
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c 



c 



1 + YVVV*V**3 + YRRR*R**3 + YDDD-'D 

1 + FY )/ (MAS-YVDOT) 

RDOT= ( NV-’-'V + NR-'R + ND*D + NVVR' ; 'V 
1 + N V RR “ V " R " “ 2 

1 + NVVV*V**3 + NRRR*R**3 + NDDD-'D* 
1 + MZ ) / (IZ-NRDOT ) 

WHEN TO PRINTOUT 

IF (ICOUNT.EQ.2 ) GO TO 50 
GO TO 300 

CONVERT RADIANS TO DEGREES 
5 0 YAWDEG = YAW*5 7.296 



*3 

*2*R 



3 



RDEG=R*57 . 296 



RDDEG=RDOT*57 .296 



DDEG=D*57.296 
YAWC=YAWC*57 .296 
WRITE (6,100) TIME, YAWDEG 
100 FORMAT ( IX , F 12 . 8 , IX , F 12 . 8 ) 
ICOUNT= 1 

C TEST IF WANT TO STOP 
300 IF (TIME . GE . ETIME ) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELT= 1 . 

C INTEGRATION 



U=U+UDOT*DELT 

V=V+VDOT*DELT 

R=R+RDOT' v DELT 

YAW=YAW+R-DELT 



X2=X2 + DX2“'DELT 
X3=X3+DX3*DELT 



C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT=U*DCOS ( YAW) - V*DSIN ( YAW) 
YDOT=U*DSIN(YAW)+V*DCOS(YAW) 
X=X+XDOT*DELT 
Y=Y+YDOT*DELT 



TIME = TIME + DELT 
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IC0UNT=IC0UNT>1 
ISE=ISE + LAMDA*YAWE**2 
ISR=ISR + D**2 
GO TO 200 

C J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 

WRITE( 6 , 500 ) ISE , I SR , TDIFF 
500 FORMAT ( ' 1 ' , 5X , ' ISE= ' , F15 . 7 , ' ISR= ' , F15 . 7 , ’ 
1 TOTAL= ' , F15 . 7 ) 

STOP 

END 

//GO.SYSIN DD * 

/* 

/ /GO . FT12F001 DD DISP=SHR,DSN=MSS . S2160 . A211 
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